Implementing the Rasch Model in Python: Fix the Scale Before Fitting

An optimizer can return plausible ability and difficulty values while the model still has no unique origin. In a Rasch model, adding the same constant to every person's ability and every item's difficulty leaves every response probability unchanged. More iterations cannot tell the optimizer which origin you intended.

This note builds a small, penalized Rasch estimator with NumPy and SciPy. The experiment uses an explicit zero-mean constraint on item difficulty, checks the gradient numerically, and includes people who answer everything correctly or incorrectly. It is a study of estimation mechanics, not a validated scoring system for students.

The probability is about a difference

For person i and item j, the model is P(correct) = sigmoid(theta[i] - b[j]). Ability is theta; difficulty is b. At equal ability and difficulty, the probability is one half. Raising ability increases the probability; raising difficulty decreases it. All items share the same discrimination, fixed at one in this parameterization.

The Stan guide's Rasch discussion provides the statistical context. The Python implementation here is a penalized joint fit, not Stan's Bayesian example and not conditional maximum likelihood. That distinction matters when interpreting the returned numbers.

Consider ability 2 and difficulty 1. The logit is 1. Ability 9 and difficulty 8 give the same logit. The likelihood cannot prefer either pair. The test suite adds seven to both fitted arrays and verifies unchanged probabilities. A successful probability calculation therefore does not establish identified parameters.

Choose an origin, then state the penalty

The implementation optimizes all person abilities and only the first J - 1 item difficulties. It derives the final difficulty as the negative sum of the others. Consequently, the item difficulties sum to zero throughout optimization; this is not a cosmetic subtraction performed after fitting.

There is a second problem. A person with all correct responses can keep increasing their unpenalized likelihood by sending ability toward positive infinity. All incorrect responses cause the opposite behavior. Our fixture deliberately contains both cases. A finite optimizer stopping point would not make those unpenalized estimates finite maximum-likelihood solutions.

The example instead adds 0.5 * penalty * (sum(theta²) + sum(b²)), with a positive finite penalty. This shrinks both parameter groups and produces finite estimates for the fixture. It changes the estimator. A stronger penalty moves the abilities toward zero, which the tests also check. The chosen value of one is an experiment setting, not a calibrated recommendation for an assessment.

Here is the objective and its analytic gradient from the executable file:

def objective(parameters, y, observed, penalty):
    theta, b = unpack(parameters, *y.shape)
    logits = theta[:, None] - b[None, :]
    residual = np.where(observed, expit(logits) - np.nan_to_num(y), 0.)
    loss = np.logaddexp(0., logits[observed]).sum()
    loss -= (y[observed] * logits[observed]).sum()
    loss += 0.5 * penalty * (theta @ theta + b @ b)
    theta_gradient = residual.sum(axis=1) + penalty * theta
    b_gradient = -residual.sum(axis=0) + penalty * b
    gradient = np.r_[theta_gradient, b_gradient[:-1] - b_gradient[-1]]
    return float(loss), gradient

The loss uses logaddexp(0, logit) - response * logit, avoiding a separate logarithm of a probability rounded to zero. The last gradient expression accounts for the dependent final difficulty. Forgetting that subtraction gives an optimizer a gradient for a different problem. A numerical gradient check is particularly useful here because a fitting curve can look reasonable despite that mistake.

A deliberately small response matrix

The fixture has five people and four items. Its rows contain four, three, two, one, and zero correct answers. Correct answers start at the first item, so the item totals also impose a clear expected difficulty order. This construction makes a sign error easy to detect; it is not a realistic sample of response behavior.

With Python 3.11.5, NumPy 2.2.6 and SciPy 1.15.3, the recorded fit returns abilities [1.0719, 0.5220, 0.0000, -0.5220, -1.0719] and difficulties [-0.7260, -0.2388, 0.2388, 0.7260], rounded to four decimals. The mean difficulty is zero. Those values belong to this matrix, penalty and constraint; they are not an external proficiency scale.

Download the implementation, tests and pinned requirements. The reproduction instructions include environment setup. After installing the requirements, run:

python3 -B -m unittest -v test_example.py

Seven tests pass. Alongside the expected ordering and gradient check, they cover missing observations, invalid responses, nonpositive penalties and forced optimizer nonconvergence. The code checks SciPy's success flag and finite parameters before returning a fit; see the versioned minimize contract.

Missing is not incorrect

A NaN marks an unobserved response and contributes no likelihood term. Replacing it with zero would assert an incorrect answer. Every person and item must have at least one observation; an entirely empty row or column is rejected rather than quietly estimated from the penalty alone.

That input rule is still weaker than a useful assessment design. Disconnected groups of people and items can obtain finite penalized estimates without observations linking their scales. Regularization supplies information that the responses do not. A practical dataset needs a connected design, uncertainty assessment and checks that item behavior fits the intended construct.

Likewise, the experiment does not test local independence, common discrimination or stability across groups. It supplies no confidence intervals and no evidence for an adaptive-testing stopping rule. Before using a Rasch fit to compare people, inspect those assumptions and the source of the shared scale. The useful engineering boundary is concrete: an identified, numerically checked objective is necessary for a trustworthy implementation, but it cannot validate the measurement problem on its own.