create_synthetic_irt_polytomous

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