Repository: eribean/girth_mcmc (language: python)

Issue:
Add support for missing values.
Need to support datasets that have missing values.


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml
new file mode 100644
index 0000000..ee222d0
--- /dev/null
+++ b/.github/workflows/python-publish.yml
@@ -0,0 +1,28 @@
+name: Girth MCMC Upload
+
+on:
+  release:
+    types: [published]
+
+jobs:
+  deploy:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: '3.x'
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install build
+    - name: Build package
+      run: python -m build
+    - name: Publish package
+      uses: pypa/gh-action-pypi-publish@27b31702a0e7fc50959f5ad993c78deac1bdfc29
+      with:
+        user: __token__
+        password: ${{ secrets.PYPI_API_TOKEN }}
\ No newline at end of file
diff --git a/README.md b/README.md
index 22bd4bf..87b2df4 100644
--- a/README.md
+++ b/README.md
@@ -13,8 +13,20 @@ packages can be installed through pip otherwise.
 * Python >= 3.7.6
 * Numpy
 * Scipy
+* Girth
 * PyMC3
 
+## Installation
+Via pip
+```
+pip install girth_mcmc --upgrade
+```
+
+From Source
+```
+pip install . -t $PYTHONPATH --upgrade
+```
+
 # Supports
 **Unidimensional**
 * Rasch Model 
@@ -27,8 +39,8 @@ packages can be installed through pip otherwise.
 Subject to change but for now:
 ```python
 import numpy as np
-from girth_mcmc import (create_synthetic_irt_dichotomous, 
-                        GirthMCMC)
+from girth import create_synthetic_irt_dichotomous
+from girth_mcmc import GirthMCMC
                         
 discrimination = 0.89 * np.sqrt(-2 * np.log(np.random.rand(10)))
 difficulty = np.random.randn(10)
@@ -46,8 +58,8 @@ print(results)
 for the graded response model, pass in the number of categories
 ```python
 import numpy as np
-from girth_mcmc import (create_synthetic_irt_polytomous, 
-                        GirthMCMC)
+from girth import create_synthetic_irt_polytomous
+from girth_mcmc import GirthMCMC
 
 n_categories = 3
 
@@ -65,13 +77,36 @@ results = girth_model(syn_data)
 print(results)
 ```
 
+Is some data missing? Tag it with a convenience function and run it like normal
+```python
+import numpy as np
+from girth import create_synthetic_irt_dichotomous
+from girth_mcmc import GirthMCMC
+from girth_mcmc.utils import tag_missing_data_mcmc
+                        
+discrimination = 0.89 * np.sqrt(-2 * np.log(np.random.rand(10)))
+difficulty = np.random.randn(10)
+theta = np.random.randn(100)
+
+syn_data = create_synthetic_irt_dichotomous(difficulty, discrimination, 
+                                            theta)
+mask = np.random.rand(*syn_data.shape) < .1
+syn_data[mask] = -9999
+syn_data_missing = tag_missing_data_mcmc(syn_data, [0, 1])
+
+girth_model = GirthMCMC(model='2PL', 
+                        options={'n_processors': 4})
+results = girth_model(syn_data_missing)
+print(results)
+```
+
 Don't like waiting? me either. Run Variational Inference for faster
 but less accurate estimation.
 
 ```python
 import numpy as np
-from girth_mcmc import (create_synthetic_irt_polytomous, 
-                        GirthMCMC)
+from girth import create_synthetic_irt_polytomous
+from girth_mcmc import GirthMCMC
 
 n_categories = 3
 
@@ -92,18 +127,47 @@ print(results_variational)
 ```
 
 ## Unittests
-The unittests are just smoke tests for now:
 
-**Without** coverage.py module
+**pytest** with coverage.py module
 ```
-nosetests testing/
+pytest --cov=girth_mcmc --cov-report term
 ```
 
-**With** coverage.py module
+**nose** with coverage.py module
 ```
-nosetests --with-coverage --cover-package=girth_mcmc testing/
+nosetests --with-coverage --cover-package=girth_mcmc
 ```
 
+## Contact
+
+Ryan Sanchez  
+ryan.sanchez@gofactr.com
+
+
 ## Other Estimation Packages
 If you are looking for Marginal Maximum Likelihood estimation routines, 
 check out [GIRTH](https://eribean.github.io/girth/).
+
+## License
+
+MIT License
+
+Copyright (c) 2021 Ryan C. Sanchez
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
\ No newline at end of file
diff --git a/girth_mcmc/__init__.py b/girth_mcmc/__init__.py
index bd33dfe..14a0f29 100644
--- a/girth_mcmc/__init__.py
+++ b/girth_mcmc/__init__.py
@@ -1,6 +1,1 @@
-from .synthetic import (create_correlated_abilities,
-                        create_synthetic_irt_dichotomous,
-                        create_synthetic_mirt_dichotomous,
-                        create_synthetic_irt_polytomous)
-
 from .girth_class import GirthMCMC
\ No newline at end of file
diff --git a/girth_mcmc/synthetic.py b/girth_mcmc/synthetic.py
deleted file mode 100644
index 080f1a9..0000000
--- a/girth_mcmc/synthetic.py
+++ /dev/null
@@ -1,300 +0,0 @@
-import numpy as np
-from scipy.special import expit
-
-
-def irt_evaluation(difficulty, discrimination, thetas):
-    """ Evaluation of unidimensional IRT model.
-
-    Evaluates an IRT model and returns the exact values.  This function
-    supports only unidimemsional models
-
-    Assumes the model
-        P(theta) = 1.0 / (1 + exp(discrimination * (theta - difficulty)))
-
-    Args:
-        difficulty: (1d array) item difficulty parameters
-        discrimination:  (1d array | number) item discrimination parameters
-        thetas: (1d array) person abilities
-
-    Returns:
-        probabilities: (2d array) evaluation of sigmoid for given inputs
-    """
-    # If discrimination is a scalar, make it an array
-    if np.atleast_1d(discrimination).size == 1:
-        discrimination = np.full_like(difficulty, discrimination,
-                                      dtype='float')
-
-    kernel = thetas - difficulty[:, None]
-    kernel *= discrimination[:, None]
-    return expit(kernel)
-
-
-def create_correlated_abilities(correlation_matrix, n_participants):
-    """ Creates correlated ability parameters based on an input correlation matrix.
-
-    This is a helper function for use in synthesizing multi-dimensional data
-    assuming multivariate normal distribution
-
-    Args:
-        correlation_matrix: (2d array) Symmetric matrix defining
-                            the correlation between the abilities
-        n_participants: number of participants to synthesize
-
-    Returns:
-        abilities: (2d array) correlated abilities
-    """
-    lower = np.linalg.cholesky(correlation_matrix)
-
-    return lower @ np.random.randn(correlation_matrix.shape[0], n_participants)
-
-
-def create_synthetic_irt_dichotomous(difficulty, discrimination, thetas,
-                                     guessing=0, seed=None):
-    """ Creates dichotomous unidimensional synthetic IRT data.
-
-    Creates synthetic IRT data to test parameters estimation functions.  
-    Only for use with dichotomous outputs
-
-    Assumes the model
-        P(theta) = 1.0 / (1 + exp(discrimination * (theta - difficulty)))
-
-    Args:
-        difficulty: [array] of difficulty parameters
-        discrimination:  [array | number] of discrimination parameters
-        thetas: [array] of person abilities
-        guessing: [array | number] of guessing parameters associated with items
-        seed: Optional setting to reproduce results
-
-    Returns:
-        synthetic_data: (2d array) realization of possible response given parameters
-
-    """
-    if seed:
-        np.random.seed(seed)
-
-    if np.ndim(guessing) < 1:
-        guessing = np.full_like(difficulty, guessing)
-
-    continuous_output = irt_evaluation(difficulty, discrimination, thetas)
-
-    # Add guessing parameters
-    continuous_output *= (1.0 - guessing[:, None])
-    continuous_output += guessing[:, None]
-
-    # convert to binary based on probability
-    random_compare = np.random.rand(*continuous_output.shape)
-
-    return (random_compare <= continuous_output).astype('int')
-
-
-def create_synthetic_mirt_dichotomous(difficulty, discrimination, thetas,
-                                      seed=None):
-    """ Creates dichotomous multidimensional synthetic IRT data.
-
-    Assumes the model
-        P(theta) = 1.0 / (1 + exp(-1 * (dot(discrimination,theta) + difficulty)))
-
-    Args:
-        difficulty: [array, M] of difficulty parameters
-        discrimination:  [2-D array, MxN] of discrimination parameters
-        thetas: [2-D array, NxP] of person abilities
-        seed: Optional setting to reproduce results
-
-    Returns:
-        synthetic_data: (2d array) realization of possible response given parameters
-
-
-    Example:
-        n_factors = 3
-        n_items = 15
-        n_people = 500
-        difficulty = np.linspace(-2.5, 2.5, n_items)
-        discrimination = np.random.randn(n_items, n_factors)
-        thetas = np.random.randn(n_factors, n_people)
-
-        synthetic_data = create_synthetic_mirt_dichotomous(difficulty, discrimination, thetas)
-    """
-    if seed:
-        np.random.seed(seed)
-
-    # If the input is just a vector of discriminations
-    if (np.ndim(discrimination) == 1) or (discrimination.shape[0] == 1):
-        discrimination = np.vstack((discrimination,) * difficulty.shape[0])
-
-    # Inline computation of the logistic kernel
-    kernel_terms = discrimination @ thetas
-    kernel_terms += difficulty[:, None]
-    continuous_output = expit(kernel_terms)
-
-    # convert to binary based on probability
-    random_compare = np.random.rand(*continuous_output.shape)
-
-    return (random_compare <= continuous_output).astype('int')
-
-
-# Private functions for polytomous outputs
-def _my_digitize(the_input):
-    """
-    Private function to compute polytomous levels.
-    The input has been concatenated to use the
-    vectorize functions (value, thresholds)       
-    """
-    return np.searchsorted(the_input[1:], the_input[0])
-
-
-def _graded_func(difficulty, discrimination, thetas, output):
-    """
-    Private function to compute the probabilities for
-    the graded response model.  This is done in place
-    and does not return anything
-    """
-    # This model is based on the difference of standard
-    # logistic functions.
-
-    # Do first level
-    output[0] = 1.0 - irt_evaluation(np.array([difficulty[0]]),
-                                     discrimination, thetas)
-
-    for level_ndx in range(1, output.shape[0]-1):
-        right = irt_evaluation(np.array([difficulty[level_ndx]]),
-                               discrimination, thetas)
-        left = irt_evaluation(np.array([difficulty[level_ndx-1]]),
-                              discrimination, thetas)
-        output[level_ndx] = left - right
-
-    # Do last level
-    output[-1] = irt_evaluation(np.array([difficulty[-1]]),
-                                discrimination, thetas)
-
-
-def _credit_func(difficulty, discrimination, thetas, output):
-    """
-    Private function to compute the probabilities for
-    the partial credit model.  This is done in place
-    and does not return anything
-    """
-    # This model is based on exponentials and normalized to
-    # make sure the expected probablity is equal to one
-    output *= 0.0  # clear any previous values
-    output[1:, :] += thetas
-    output[1:, :] -= difficulty[:, None]
-    output *= discrimination
-    np.cumsum(output, axis=0, out=output)
-    np.exp(output, out=output)
-
-    normalizing_term = 1.0 / np.sum(output, axis=0)
-    output *= normalizing_term
-
-
-def _unfold_func(difficulty, discrimination, thetas, output, src=0.):
-    """
-    Private function to compute the probabilities for
-    the graded unfolding model.  This is done in place
-    and does not return anything
-    """
-    # Call partial credit model on difficulties
-    _credit_func(difficulty, discrimination, thetas, output)
-
-    # Add the probabilities together, should always be even
-    for ndx in range(output.shape[0] // 2):
-        output[ndx] += (output[-(ndx+1)] * (1. - src))
-        output[-(ndx+1)] *= src
-
-
-def _check_difficulty_parameters(difficulty, model):
-    """
-    Given a model type, check the difficulty parameters
-    for validity
-    """
-    max_value = difficulty.shape[1] + 1
-
-    if model in ["grm"]:
-        # Check that all the arguments are sorted
-        if not np.all(difficulty[:, :-1] < difficulty[:, 1:]):
-            raise AssertionError("Difficulty Parameters must be "
-                                 "in ascending order")
-
-    elif model in ['gum']:
-        # Parameters must be odd
-        if max_value % 2:
-            raise AssertionError("There must be an odd number of "
-                                 "difficulty parameters")
-
-        # Parameters must be skew-symmetric about the center point
-        middle_index = (difficulty.shape[1] - 1) // 2
-        adjusted_difficulty = (difficulty -
-                               difficulty[:, middle_index][:, None])
-
-        if not np.all(np.abs(adjusted_difficulty.sum(axis=1)) < 1e-7):
-            raise AssertionError("Difficulty Parameters must be "
-                                 "symmetric about offset")
-
-        max_value = middle_index + 1
-
-    return max_value
-
-
-def create_synthetic_irt_polytomous(difficulty, discrimination, thetas,
-                                    model='grm', seed=None):
-    """ Creates polytomous unidimensional synthetic IRT data.
-
-    Creates polytomous output with specified number of levels from [1, levels]
-
-    Args:
-        difficulty: [2D array (items x n_levels-1)] of difficulty parameters
-        discrimination:  [array | number] of discrimination parameters
-        thetas: [array] of person abilities
-        model: ["grm", "pcm", "gum"] string specifying which polytomous model to use
-                'grm': Graded Response Model
-                'pcm': Generalized Partial Credit Model
-                'gum': Generalized Graded Unfolding Model
-        seed: Optional setting to reproduce results
-
-    Returns:
-        synthetic_data: (2d array) realization of possible response given parameters
-    """
-    difficulty = np.atleast_2d(difficulty)
-    n_items, n_levels = difficulty.shape
-
-    if n_levels == 1:
-        raise AssertionError("Polytomous items must have more than 1 threshold")
-
-    if seed:
-        np.random.seed(seed)
-
-    # Check for single input of discrimination
-    if np.atleast_1d(discrimination).size == 1:
-        discrimination = np.full((n_items,), discrimination)
-
-    # Get the model to use, will throw error if not supported
-    probability_func = {'grm': _graded_func,
-                        'pcm': _credit_func,
-                        'gum': _unfold_func}[model.lower()]
-
-    # Check difficulty parameters for validity
-    clip_high = _check_difficulty_parameters(difficulty, model.lower())
-
-    # Initialize output for memory concerns
-    level_scratch = np.zeros((n_levels + 2, thetas.size))
-    output = np.zeros((n_items, thetas.size), dtype='int')
-
-    # Loop over items and compute probability estimates
-    # for each of the levels and assign level based on
-    # those probabilities
-    for item_ndx in range(n_items):
-        # Obtain the probabilities for the data (in-place)
-        probability_func(difficulty[item_ndx], discrimination[item_ndx],
-                         thetas, level_scratch[1:, :])
-
-        # Get the thresholds of the levels
-        np.cumsum(level_scratch[1:, :], axis=0, out=level_scratch[1:, :])
-        level_scratch[0] = np.random.rand(thetas.size)
-
-        # Discritize the outputs based on the thresholds
-        output[item_ndx] = np.apply_along_axis(
-            _my_digitize, axis=0, arr=level_scratch)
-
-    # Add 1 to return [1, n_levels]
-    output += 1
-    np.clip(output, 1, clip_high, out=output)
-    return output
diff --git a/girth_mcmc/utils/__init__.py b/girth_mcmc/utils/__init__.py
index 389934d..69aae8d 100644
--- a/girth_mcmc/utils/__init__.py
+++ b/girth_mcmc/utils/__init__.py
@@ -1,2 +1,3 @@
 from .options import validate_mcmc_options
-from .rayleigh import Rayleigh
\ No newline at end of file
+from .missing_data import *
+from .rayleigh import Rayleigh
diff --git a/girth_mcmc/utils/missing_data.py b/girth_mcmc/utils/missing_data.py
new file mode 100644
index 0000000..bb65461
--- /dev/null
+++ b/girth_mcmc/utils/missing_data.py
@@ -0,0 +1,22 @@
+from numpy import ma, isin
+
+
+__all__ = ['tag_missing_data_mcmc']
+
+
+def tag_missing_data_mcmc(dataset, valid_responses):
+    """Checks the data for valid responses.
+    
+    Args:
+        dataset: (array) array to validate
+        valid_responses: (array-like) list of valid responses
+        
+    Returns:
+        updated_dataset: (array) data that holds only valid_responses and
+                         invalid_fill
+    """
+    mask = isin(dataset, valid_responses)
+
+    # MCMC uses a masked array to identify missing data
+    return ma.masked_array(dataset, ~mask)
+
diff --git a/requirements.txt b/requirements.txt
index c398808..e840506 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,4 @@
 numpy>=1.18.2
 scipy>=1.4.1
 pymc3>=3.10.0
+girth>=0.5.0
diff --git a/setup.py b/setup.py
index 15087a8..7563bd4 100644
--- a/setup.py
+++ b/setup.py
@@ -1,30 +1,27 @@
 from setuptools import setup, convert_path
-
-# read the contents of your README file
 from os import path
+
 this_directory = path.abspath(path.dirname(__file__))
 with open(path.join(this_directory, 'README.md'), encoding='utf-8') as f:
     long_description = f.read()
 
 
-print(convert_path('/performance'))
 if __name__ == '__main__':
     setup(
         name="girth_mcmc", 
         packages=['girth_mcmc', 'girth_mcmc.dichotomous', 'girth_mcmc.polytomous', 
                   'girth_mcmc.utils'],
         package_dir={'girth_mcmc': 'girth_mcmc'},
-        version="0.1.1",
+        version="0.2.0",
         license="MIT",
-        description="MCMC Item Response Theory Estimation.",
+        description="Bayesian Item Response Theory Estimation.",
         long_description=long_description.replace('<ins>','').replace('</ins>',''),
         long_description_content_type='text/markdown',
         author='Ryan C. Sanchez',
         author_email='rsanchez44@gatech.edu',
         url = 'https://github.com/eribean/girth_mcmc/',
-        download_url = 'https://github.com/eribean/girth_mcmc/archive/v0.1.1.tar.gz',
-        keywords = ['IRT', 'Psychometrics', 'Item Response Theory', 'MCMC'],
-        install_requires = ['numpy', 'scipy', 'pymc3'],
+        keywords = ['IRT', 'Psychometrics', 'Item Response Theory', 'MCMC', 'Bayesian'],
+        install_requires = ['numpy', 'scipy', 'pymc3', 'girth'],
         classifiers = [
             'Development Status :: 3 - Alpha',
             'Intended Audience :: Science/Research',
@@ -32,5 +29,6 @@
             'License :: OSI Approved :: MIT License',
             'Programming Language :: Python :: 3.7',
             'Programming Language :: Python :: 3.8'
+            'Programming Language :: Python :: 3.9'
         ]
     )


Test patch:
diff --git a/testing/test_dichotomous.py b/test/test_dichotomous.py
similarity index 98%
rename from testing/test_dichotomous.py
rename to test/test_dichotomous.py
index 1355375..515c946 100644
--- a/testing/test_dichotomous.py
+++ b/test/test_dichotomous.py
@@ -2,9 +2,8 @@
 
 import numpy as np
 
-from girth_mcmc import (create_synthetic_irt_dichotomous, 
-                        GirthMCMC)
-
+from girth import create_synthetic_irt_dichotomous
+from girth_mcmc import GirthMCMC
 
 
 class TestDichotomous(unittest.TestCase):
diff --git a/testing/test_polytomous.py b/test/test_polytomous.py
similarity index 95%
rename from testing/test_polytomous.py
rename to test/test_polytomous.py
index 0453feb..479aec2 100644
--- a/testing/test_polytomous.py
+++ b/test/test_polytomous.py
@@ -2,8 +2,8 @@
 
 import numpy as np
 
-from girth_mcmc import (create_synthetic_irt_polytomous, 
-                        GirthMCMC)
+from girth import create_synthetic_irt_polytomous
+from girth_mcmc import GirthMCMC
 
 
 class TestPolytomous(unittest.TestCase):
diff --git a/testing/test_synthetic.py b/testing/test_synthetic.py
deleted file mode 100644
index fce3535..0000000
--- a/testing/test_synthetic.py
+++ /dev/null
@@ -1,337 +0,0 @@
-import unittest
-
-import numpy as np
-
-from girth_mcmc import create_synthetic_irt_dichotomous
-from girth_mcmc import create_synthetic_mirt_dichotomous
-from girth_mcmc import create_synthetic_irt_polytomous
-from girth_mcmc import create_correlated_abilities
-
-from girth_mcmc.synthetic import (_my_digitize, _credit_func, 
-                                  _graded_func, _unfold_func,
-                                  _check_difficulty_parameters)
-
-
-class TestSynthetic(unittest.TestCase):
-    """Testing the creation of synthetic irt function."""
-
-    def test_synthetic_irt_creation(self):
-        """Testing the creation of synthetic data."""
-        seed = 31
-
-        # Regression test
-        expected = np.array([[False, False, False, False,  True,  True],
-                             [False, False,  True,  True,  True,  True],
-                             [False, False, False,  True,  True,  True]])
-
-        value = create_synthetic_irt_dichotomous(np.array([1.2, -0.2, 1.3]),
-                                                 1.31, np.linspace(-6, 6, 6),
-                                                 seed=seed)
-
-        np.testing.assert_array_equal(expected, value)
-
-
-    def test_synthetic_mirt_creation(self):
-        """Testing the creation of synthetic data."""
-        seed = 164
-        np.random.seed(seed-1)
-        # Regression test
-        expected = np.array([[False, False, False, False, False, False],
-                             [ True, False, False,  True, False,  True],
-                             [ True,  True, False, False, False, False],
-                             [ True,  True,  True,  True, False,  True],
-                             [ True,  True,  True,  True,  True,  True]])
-
-        n_factors = 3
-        n_items = 5
-        n_people = 6
-        discrimination = np.random.randn(n_items, n_factors)
-        difficulty = np.linspace(-5, 5, n_items)
-        thetas = np.random.randn(n_factors, n_people)
-        value = create_synthetic_mirt_dichotomous(difficulty, discrimination,
-                                                  thetas, seed)
-
-        np.testing.assert_array_equal(expected, value)
-
-
-    def test_synthetic_mirt_creation_single(self):
-        """Testing the creation of synthetic data, common discrimination."""
-        seed = 546
-        np.random.seed(seed-1)
-        # Regression test
-        expected = np.array([[False, False, False, False, False, False],
-                             [False, False, False, False, False, False],
-                             [False, False, False,  True,  True,  True],
-                             [False,  True,  True,  True,  True,  True],
-                             [ True,  True,  True,  True,  True,  True]])
-
-        n_factors = 3
-        n_items = 5
-        n_people = 6
-        discrimination = np.random.randn(1, n_factors)
-        difficulty = np.linspace(-5, 5, n_items)
-        thetas = np.random.randn(n_factors, n_people)
-        value = create_synthetic_mirt_dichotomous(difficulty, discrimination,
-                                                  thetas, seed)
-
-        np.testing.assert_array_equal(expected, value)
-
-
-    def test_correlated_abilities(self):
-        """Testing the creation of correlated abilities."""
-        np.random.seed(120)
-        n_participants = 1000
-        rho = 0.73
-        correlation_matrix = np.array([[1, rho], [rho, 1]])
-
-        output = create_correlated_abilities(correlation_matrix, n_participants)
-        output_corr = np.corrcoef(output)
-
-        np.testing.assert_almost_equal(output_corr, correlation_matrix, decimal=1)
-
-
-class TestPolytomousSynthetic(unittest.TestCase):
-    """Testing the creation of synthetic polytomous irt function."""
-
-    def test_my_digititize(self):
-        """Testing local digitize function."""
-        test_candidate = np.zeros((3,))
-
-        # Zero position
-        test_candidate[:] = [0.25, .3, .6]
-        output = _my_digitize(test_candidate)
-
-        self.assertEqual(output, 0)
-
-        # One position
-        test_candidate[:] = [0.3, .25, .6]
-        output = _my_digitize(test_candidate)
-
-        self.assertEqual(output, 1)
-
-        # One position boundary
-        test_candidate[:] = [0.3, .25, .3]
-        output = _my_digitize(test_candidate)
-
-        self.assertEqual(output, 1)
-
-        # Two position 
-        test_candidate[:] = [0.35, .25, .3]
-        output = _my_digitize(test_candidate)
-
-        self.assertEqual(output, 2)
-
-
-    def test_check_difficulty_parameters(self):
-        """Testing the validity check of polytomous difficulty."""
-        np.random.seed(513)
-        n_size = 5
-        difficulty = np.random.rand(5, n_size)
-
-        # PCM
-        result = _check_difficulty_parameters(np.sort(difficulty, axis=1), 'grm')
-        self.assertEqual(result, n_size + 1)
-
-        # GRM
-        with self.assertRaises(AssertionError):
-            _check_difficulty_parameters(difficulty, 'grm')
-
-        result = _check_difficulty_parameters(np.sort(difficulty, axis=1), 'grm')
-        self.assertEqual(result, n_size + 1)
-
-        #GUM
-        difficulty = np.arange(n_size*n_size).reshape(n_size, n_size)
-
-        result = _check_difficulty_parameters(difficulty, 'gum')
-        self.assertEqual(result, n_size // 2 + 1)
-
-        with self.assertRaises(AssertionError):
-            _check_difficulty_parameters(difficulty[:, :-1], 'gum')
-
-        difficulty[:, 0] = np.random.permutation(difficulty[:, 0])
-        with self.assertRaises(AssertionError):
-            _check_difficulty_parameters(difficulty, 'gum')
-
-
-    def test_graded_function(self):
-        """Testing the graded response model computation"""
-        # Create basic data
-        difficulties = np.array([-2.3, .3, 1.2])
-        discrimination = 0.78
-        thetas = np.linspace(-3, 3, 100)
-
-        # Initialize output variable and call data
-        output = np.zeros((difficulties.size + 1, thetas.size))
-        _graded_func(difficulties, discrimination, thetas, output)
-
-        # Compare to hand computations
-        first_position = 1.0 / (1.0 + np.exp(discrimination * (thetas - difficulties[0])))
-        second_position = 1.0 / (1.0 + np.exp(discrimination * (thetas - difficulties[1])))
-        last_position = 1.0 / (1.0 + np.exp(discrimination * (thetas - difficulties[2])))
-
-        np.testing.assert_array_almost_equal(output[0], first_position)
-        np.testing.assert_array_almost_equal(output[1], second_position - first_position)
-        np.testing.assert_array_almost_equal(output[2], last_position - second_position)
-        np.testing.assert_array_almost_equal(output[3], 1 - last_position)
-        
-
-    def test_credit_function(self):
-        """Testing the partial credit computation"""
-        # Create basic data
-        difficulties = np.array([-0.67, .24, .84])
-        discrimination = 1.24
-        thetas = np.linspace(-3, 3, 100)
-
-        # Initialize output variable and call data
-        output = np.zeros((difficulties.size + 1, thetas.size))
-        _credit_func(difficulties, discrimination, thetas, output)
-
-        # Compare to hand computations
-        first_position = 1.0 
-        second_position = np.exp(discrimination * (thetas - difficulties[0]))
-        third_position = np.exp(discrimination * (thetas - difficulties[1])) * second_position
-        last_position = np.exp(discrimination * (thetas - difficulties[2])) * third_position
-
-        normalizing = first_position + second_position + third_position + last_position
-
-        np.testing.assert_array_almost_equal(output[0], first_position / normalizing)
-        np.testing.assert_array_almost_equal(output[1], second_position / normalizing)
-        np.testing.assert_array_almost_equal(output[2], third_position / normalizing)
-        np.testing.assert_array_almost_equal(output[3], last_position / normalizing)
-
-
-    def test_unfolding_function(self):
-        """Testing the graded unfolding model computation."""
-        np.random.seed(831)
-        offset = 0.3
-        difficulties = np.array([-1.3, -.7, -.3, 0, .3, .7, 1.3])
-        discrimination = 0.871
-        thetas = np.linspace(-3, 3, 100)
-        
-        # Initialize output variable and call data
-        output = np.zeros((difficulties.size + 1, thetas.size))
-        _unfold_func(difficulties + offset, discrimination, thetas, output)
-
-        # Compare to hand calculations
-        span = difficulties.size / 2
-        beta = discrimination * (thetas - offset)
-        first_position = (np.cosh((span - 0) * beta) * 
-                          np.exp(-discrimination * 0))
-        second_position = (np.cosh((span - 1) * beta) * 
-                           np.exp(-discrimination * difficulties[:1].sum()))
-        third_position = (np.cosh((span - 2) * beta) * 
-                           np.exp(-discrimination * difficulties[:2].sum()))
-        last_position = (np.cosh((span - 3) * beta) * 
-                         np.exp(-discrimination * difficulties[:3].sum()))
-
-        normalizing = first_position + second_position + third_position + last_position
-
-        np.testing.assert_array_almost_equal(output[0], first_position / normalizing)
-        np.testing.assert_array_almost_equal(output[1], second_position / normalizing)
-        np.testing.assert_array_almost_equal(output[2], third_position / normalizing)
-        np.testing.assert_array_almost_equal(output[3], last_position / normalizing)
-
-        # Test SRC works too
-        _unfold_func(difficulties + offset, discrimination, thetas, output, src=1.)
-        np.testing.assert_array_almost_equal(output[0] + output[-1], first_position / normalizing)
-        np.testing.assert_array_almost_equal(output[1] + output[-2], second_position / normalizing)
-        np.testing.assert_array_almost_equal(output[2] + output[-3], third_position / normalizing)
-        np.testing.assert_array_almost_equal(output[3] + output[-4], last_position / normalizing)
-
-
-    def test_create_polytomous_data_fail(self):
-        """Testing synthetic polytomous function fails with 1 level"""
-        np.random.seed(808)
-        difficulty = np.array([[1.0]])
-        bad_graded_difficulty = np.random.rand(4, 5)
-        even_gum_difficulty = np.random.rand(4, 4)
-
-        with self.assertRaises(AssertionError):
-            create_synthetic_irt_polytomous(difficulty, difficulty, difficulty)
-
-        with self.assertRaises(KeyError):
-            create_synthetic_irt_polytomous([1, 2, 3], difficulty, 
-                                            difficulty, model='boom')
-
-        with self.assertRaises(AssertionError):
-            create_synthetic_irt_polytomous(bad_graded_difficulty, 1.2,
-                                            np.random.randn(100))
-        
-        with self.assertRaises(AssertionError):
-            create_synthetic_irt_polytomous(even_gum_difficulty, 1.2,
-                                            np.random.randn(100),
-                                            model="gum")
-
-        # Test for skew symmetric fail
-        with self.assertRaises(AssertionError):
-            create_synthetic_irt_polytomous(bad_graded_difficulty, 1.2,
-                                            np.random.randn(100),
-                                            model="gum")
-
-
-    def test_check_polytomous_discrimination(self):
-        """Smoke tests if a single value for discrimination passes"""
-        difficulty = np.array([[1.0, 2., 3.]])
-        discrimination = 3
-
-        # Simple Smoke tests
-        create_synthetic_irt_polytomous(difficulty, discrimination, difficulty)
-        create_synthetic_irt_polytomous(difficulty, np.array([discrimination]), 
-                                        difficulty)
-
-
-    def test_check_polytomous_regression(self):
-        """Regression testing graded and credit polytomous functions"""
-        seed = 876
-        np.random.seed(seed)
-        difficulty_gum = np.array([np.linspace(-1.2, 1.2, 5),
-                                   np.linspace(-.8, .8, 5),
-                                   np.linspace(-1.7, 1.7, 5),
-                                   np.linspace(-1.0, 1.0, 5)])
-        difficulty_pcm = np.random.randn(5, 4)
-        difficulty = np.sort(difficulty_pcm, axis=1)
-        discrimination = 1.23
-        thetas = np.random.randn(8)
-
-        # Regression Tests
-        poly_data_graded = create_synthetic_irt_polytomous(difficulty, 
-                                                           discrimination,
-                                                           thetas,
-                                                           model='grm',
-                                                           seed=seed)
-
-        poly_data_credit = create_synthetic_irt_polytomous(difficulty_pcm, 
-                                                           discrimination,
-                                                           thetas,
-                                                           model='pcm',
-                                                           seed=seed)
-
-        poly_data_gum = create_synthetic_irt_polytomous(difficulty_gum, 
-                                                        discrimination,
-                                                        thetas,
-                                                        model='gum',
-                                                        seed=seed)
-
-        expected_graded = np.array([[5, 2, 2, 5, 2, 2, 5, 5],
-                                    [5, 5, 1, 5, 5, 1, 2, 1],
-                                    [5, 5, 1, 3, 2, 1, 5, 1],
-                                    [3, 5, 3, 5, 4, 1, 4, 1],
-                                    [5, 5, 3, 4, 3, 3, 5, 2]])
-
-        expected_partial = np.array([[5, 5, 2, 5, 2, 1, 5, 3],
-                                     [5, 5, 1, 5, 4, 1, 4, 1],
-                                     [5, 5, 1, 3, 1, 1, 5, 1],
-                                     [5, 5, 1, 5, 4, 1, 5, 1],
-                                     [5, 5, 2, 5, 4, 2, 5, 2]])
-
-        expected_unfold = np.array([[2, 1, 2, 2, 2, 2, 2, 3],
-                                    [2, 2, 1, 2, 3, 1, 1, 1],
-                                    [2, 3, 2, 2, 3, 1, 2, 1],
-                                    [1, 3, 3, 3, 3, 1, 1, 1]])
-
-        np.testing.assert_array_equal(poly_data_graded, expected_graded)
-        np.testing.assert_array_equal(poly_data_credit, expected_partial)
-        np.testing.assert_array_equal(poly_data_gum, expected_unfold)
-
-if __name__ == '__main__':
-    unittest.main()