From 629ca2dfcf680f6a190d2f4301e9a99170517407 Mon Sep 17 00:00:00 2001 From: Alex Kholodniak Date: Sun, 11 Oct 2026 06:55:15 +0300 Subject: [PATCH 1/6] fix: make joint estimation correct and honest about convergence The 3PL gradients omitted chain-rule factors, so 3PL fits optimized the wrong objective and recovered guessing parameters at chance level. Extreme response patterns have no finite estimate, yet they were returned as ordinary numbers that kept growing with max_iter, and the latent scale was left unidentified. The stopping rule could also fire after the learning rate collapsed, far from any optimum, with no way for callers to tell. Exclude extreme and unobserved persons/items (returned as nil), fix the scale after fitting, use a stable log-sigmoid likelihood, let the learning rate recover after accepted steps, and report converged/iterations with a stationarity check (gradient_tolerance). Finite-difference and recovery specs guard the math; docs no longer claim production readiness. --- CHANGELOG.md | 14 +- README.md | 30 +++- benchmarks/convergence_benchmark.rb | 88 +++------ irt_ruby.gemspec | 15 +- lib/irt_ruby/estimation.rb | 95 ++++++++++ lib/irt_ruby/model_options_validator.rb | 3 +- lib/irt_ruby/rasch_model.rb | 135 ++++++++------ lib/irt_ruby/three_parameter_model.rb | 190 ++++++++++++-------- lib/irt_ruby/two_parameter_model.rb | 149 ++++++++------- spec/irt_ruby/estimation_spec.rb | 174 ++++++++++++++++++ spec/irt_ruby/three_parameter_model_spec.rb | 22 +-- spec/irt_ruby/two_parameter_model_spec.rb | 17 +- spec/spec_helper.rb | 1 + 13 files changed, 633 insertions(+), 300 deletions(-) create mode 100644 lib/irt_ruby/estimation.rb create mode 100644 spec/irt_ruby/estimation_spec.rb diff --git a/CHANGELOG.md b/CHANGELOG.md index 86bfbff..e5706cc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,19 @@ All notable changes to this project are documented in this file. ## [Unreleased] +### Fixed +- Fixed 3PL gradients for ability, difficulty, discrimination, and guessing; they now match finite differences of the log-likelihood. +- Respondents and items with all-correct, all-incorrect, or unobserved response patterns are excluded from estimation and returned as `nil` instead of drifting with `max_iter` or keeping random initial values. +- Log-likelihoods are computed with a numerically stable log-sigmoid instead of adding `1e-15` inside `Math.log`, which flattened the objective for badly mispredicted responses. +- Models copy response data on construction and return fresh arrays from `fit`, so later mutation of either does not change model state. + +### Changed +- `fit` fixes the latent scale: Rasch abilities are centered at 0; 2PL and 3PL abilities are standardized to mean 0 and SD 1 with likelihood-preserving item transformations. +- `fit` results include `converged:` and `iterations:`. `converged` is true only at a stationary point, controlled by the new `gradient_tolerance:` option (default `1e-3`). +- Clarified the gem summary and description to describe implemented capabilities rather than production-readiness claims. + ### Added +- Added finite-difference gradient, identification, extreme-pattern, and parameter-recovery specs. - Added executable specs for README usage examples to keep the documented Rasch, 2PL, 3PL, and missing-data flows aligned with the public API. ### Changed @@ -13,7 +25,7 @@ All notable changes to this project are documented in this file. --- -## [0.3.0] - 2025-01-14 +## [0.3.0] - 2025-06-13 ### Changed - **Code Quality** diff --git a/README.md b/README.md index f84e8a4..33487d5 100644 --- a/README.md +++ b/README.md @@ -46,8 +46,12 @@ result = model.fit # Output the estimated abilities and difficulties puts "Abilities: #{result[:abilities]}" puts "Difficulties: #{result[:difficulties]}" +puts "Converged: #{result[:converged]} after #{result[:iterations]} iterations" ``` +The third respondent answered every item correctly, so their ability has no +finite estimate and is returned as `nil` (see [Estimation Details](#estimation-details)). + Response data passed to model constructors must be either a `Matrix` or an array of arrays. Each response value must be the integer `0`, the integer `1`, or `nil` for missing data; floats such as `0.0`/`1.0`, strings, booleans, and @@ -69,6 +73,10 @@ puts three_pl_result[:discriminations] puts three_pl_result[:guessings] ``` +Three respondents are far too few for 2PL or 3PL: on this matrix `fit` returns +`converged: false`. Discrimination and guessing estimates need hundreds of +respondents. + ## Handling Missing Data Real-world data often has missing responses. Each model (Rasch, 2PL, 3PL) accepts a `missing_strategy:` option to handle nil entries: @@ -102,16 +110,20 @@ This flexibility helps you handle datasets where missingness might signify a ski ### Adaptive Learning Rate & Convergence By default, each model uses a gradient ascent with: -- An adaptive learning rate (if log-likelihood decreases, it reverts the step and reduces the rate). -- Multiple convergence checks (change in log-likelihood and average parameter updates). +- An adaptive learning rate: if log-likelihood decreases, the step is reverted and the rate is multiplied by `decay_factor`; after each accepted step the rate grows by 10%. +- A stopping rule on the change in log-likelihood and the average parameter update. +- A stationarity check that decides whether the stopping point is reported as converged. You can customize: - `max_iter`: The maximum number of iterations. -- `tolerance` and `param_tolerance`: Convergence thresholds for log-likelihood change and parameter updates. +- `tolerance` and `param_tolerance`: Stopping thresholds for log-likelihood change and parameter updates. +- `gradient_tolerance` (default `1e-3`): `fit` reports `converged: true` only if, for every parameter, the gradient averaged over the responses it depends on is below this value. - `learning_rate`: Initial learning rate. - `decay_factor`: Factor by which the learning rate is reduced on a failed step. +`fit` returns `converged: false` when `max_iter` runs out or when the stopping rule fires away from a stationary point, for example because the learning rate became too small to move. Check this flag before using the estimates; a larger `learning_rate` or `max_iter` usually helps. + Example: ```ruby IrtRuby::TwoParameterModel.new( @@ -140,10 +152,18 @@ The `seed:` keyword is available for `RaschModel`, `TwoParameterModel`, and `Thr ### Parameter Clamping For 2PL and 3PL: -- **Discriminations** (`a`) are clamped between `0.01` and `5.0`. +- **Discriminations** (`a`) are clamped between `0.01` and `5.0` during optimization. - **Guessings** (`c`, 3PL only) are clamped to `[0.0, 0.35]`. -This prevents extreme or invalid parameter estimates. +These bounds are model restrictions: a 3PL item with a true lower asymptote above `0.35` cannot be represented. Because the final scale is standardized after optimization (see below), returned discriminations can fall slightly outside the optimization bounds. + +## Estimation Details + +All models estimate person and item parameters jointly by maximizing the likelihood (joint maximum likelihood, JMLE). + +- **Extreme and unobserved patterns.** A respondent who answered every observed item correctly, or every item incorrectly, has no finite maximum-likelihood ability; the same holds for items that everyone answered correctly or incorrectly, and for rows or columns with no observed responses. These persons and items are excluded from estimation and returned as `nil`. Exclusion is repeated until stable, since removing an item can make a respondent's remaining pattern extreme. Missing-data strategies are applied first, so `:treat_as_correct` and `:treat_as_incorrect` can create extreme patterns. +- **Scale identification.** The latent scale is fixed after fitting: Rasch abilities are centered to mean 0; 2PL and 3PL abilities are standardized to mean 0 and standard deviation 1, with difficulties and discriminations transformed so that response probabilities are unchanged. Only estimable respondents enter these statistics. +- **Known limitations.** JMLE item estimates are biased when tests are short, and the bias does not vanish as the number of respondents grows. 2PL and 3PL JMLE is also slow to converge with plain gradient ascent. For high-stakes calibration, compare against marginal maximum likelihood software such as R's `mirt`. ## Performance Benchmarks diff --git a/benchmarks/convergence_benchmark.rb b/benchmarks/convergence_benchmark.rb index 589629e..435151e 100755 --- a/benchmarks/convergence_benchmark.rb +++ b/benchmarks/convergence_benchmark.rb @@ -5,46 +5,6 @@ require "irt_ruby" require "benchmark" -# Enhanced model classes that track iterations and convergence -class TrackedRaschModel < IrtRuby::RaschModel - attr_reader :iterations, :final_log_likelihood, :convergence_reason - - def fit - @iterations = 0 - prev_ll = log_likelihood - @final_log_likelihood = prev_ll - @convergence_reason = :max_iterations - - @max_iter.times do - @iterations += 1 - grad_abilities, grad_difficulties = compute_gradient - - old_a, old_d = apply_gradient_update(grad_abilities, grad_difficulties) - - current_ll = log_likelihood - param_delta = average_param_update(old_a, old_d) - - if current_ll < prev_ll - @abilities = old_a - @difficulties = old_d - @learning_rate *= @decay_factor - else - ll_diff = (current_ll - prev_ll).abs - @final_log_likelihood = current_ll - - if ll_diff < @tolerance && param_delta < @param_tolerance - @convergence_reason = :tolerance_reached - break - end - - prev_ll = current_ll - end - end - - { abilities: @abilities, difficulties: @difficulties } - end -end - def generate_data(num_people, num_items, difficulty_range: (-2..2), ability_range: (-2..2)) # Generate realistic IRT data based on known parameters true_abilities = Array.new(num_people) { rand(ability_range) } @@ -84,27 +44,27 @@ def generate_data(num_people, num_items, difficulty_range: (-2..2), ability_rang times = [] iterations = [] - convergence_reasons = [] + convergences = [] 5.times do time = Benchmark.measure do - model = TrackedRaschModel.new( + model = IrtRuby::RaschModel.new( data, max_iter: 2000, tolerance: config[:tolerance], param_tolerance: config[:param_tolerance], learning_rate: 0.01 ) - model.fit - iterations << model.iterations - convergence_reasons << model.convergence_reason + result = model.fit + iterations << result[:iterations] + convergences << result[:converged] end.real times << time end avg_time = times.sum / times.size avg_iterations = iterations.sum.to_f / iterations.size - convergence_rate = convergence_reasons.count(:tolerance_reached) / 5.0 + convergence_rate = convergences.count(true) / 5.0 printf(" Time: %6.3fs Iterations: %6.1f Convergence Rate: %4.0f%%\n", avg_time, avg_iterations, convergence_rate * 100) @@ -128,27 +88,27 @@ def generate_data(num_people, num_items, difficulty_range: (-2..2), ability_rang times = [] iterations = [] - convergence_reasons = [] + convergences = [] 5.times do time = Benchmark.measure do - model = TrackedRaschModel.new( + model = IrtRuby::RaschModel.new( data, max_iter: 1000, tolerance: 1e-5, param_tolerance: 1e-5, learning_rate: config[:rate] ) - model.fit - iterations << model.iterations - convergence_reasons << model.convergence_reason + result = model.fit + iterations << result[:iterations] + convergences << result[:converged] end.real times << time end avg_time = times.sum / times.size avg_iterations = iterations.sum.to_f / iterations.size - convergence_rate = convergence_reasons.count(:tolerance_reached) / 5.0 + convergence_rate = convergences.count(true) / 5.0 printf(" Time: %6.3fs Iterations: %6.1f Convergence Rate: %4.0f%%\n", avg_time, avg_iterations, convergence_rate * 100) @@ -171,7 +131,7 @@ def generate_data(num_people, num_items, difficulty_range: (-2..2), ability_rang times = [] iterations = [] - convergence_reasons = [] + convergences = [] 3.times do dataset = generate_data( @@ -182,23 +142,23 @@ def generate_data(num_people, num_items, difficulty_range: (-2..2), ability_rang ) time = Benchmark.measure do - model = TrackedRaschModel.new( + model = IrtRuby::RaschModel.new( dataset[:data], max_iter: 1000, tolerance: 1e-5, param_tolerance: 1e-5, learning_rate: 0.01 ) - model.fit - iterations << model.iterations - convergence_reasons << model.convergence_reason + result = model.fit + iterations << result[:iterations] + convergences << result[:converged] end.real times << time end avg_time = times.sum / times.size avg_iterations = iterations.sum.to_f / iterations.size - convergence_rate = convergence_reasons.count(:tolerance_reached) / 3.0 + convergence_rate = convergences.count(true) / 3.0 printf(" Time: %6.3fs Iterations: %6.1f Convergence Rate: %4.0f%%\n", avg_time, avg_iterations, convergence_rate * 100) @@ -233,11 +193,11 @@ def generate_data(num_people, num_items, difficulty_range: (-2..2), ability_rang times = [] iterations = [] - convergence_reasons = [] + convergences = [] 3.times do time = Benchmark.measure do - model = TrackedRaschModel.new( + model = IrtRuby::RaschModel.new( data_with_missing, max_iter: 1000, tolerance: 1e-5, @@ -245,16 +205,16 @@ def generate_data(num_people, num_items, difficulty_range: (-2..2), ability_rang learning_rate: 0.01, missing_strategy: config[:strategy] ) - model.fit - iterations << model.iterations - convergence_reasons << model.convergence_reason + result = model.fit + iterations << result[:iterations] + convergences << result[:converged] end.real times << time end avg_time = times.sum / times.size avg_iterations = iterations.sum.to_f / iterations.size - convergence_rate = convergence_reasons.count(:tolerance_reached) / 3.0 + convergence_rate = convergences.count(true) / 3.0 printf(" Time: %6.3fs Iterations: %6.1f Convergence Rate: %4.0f%%\n", avg_time, avg_iterations, convergence_rate * 100) diff --git a/irt_ruby.gemspec b/irt_ruby.gemspec index f8fead6..f3f2db4 100644 --- a/irt_ruby.gemspec +++ b/irt_ruby.gemspec @@ -8,9 +8,9 @@ Gem::Specification.new do |spec| spec.authors = ["Alex Kholodniak"] spec.email = ["alexandrkholodniak@gmail.com"] - spec.summary = "Production-ready Item Response Theory (IRT) models with comprehensive performance benchmarking and adaptive optimization." + spec.summary = "Rasch, 2PL, and 3PL Item Response Theory (IRT) models estimated by joint maximum likelihood." spec.description = <<~DESC - IrtRuby is a comprehensive Ruby library for Item Response Theory (IRT) analysis,#{" "} + IrtRuby is a Ruby library for Item Response Theory (IRT) analysis,#{" "} commonly used in educational assessment, psychological testing, and survey research. Features three core IRT models: @@ -19,13 +19,12 @@ Gem::Specification.new do |spec| • Three-Parameter Model (3PL) - Includes guessing parameter Key capabilities: - • Robust gradient ascent optimization with adaptive learning rates - • Flexible missing data strategies (ignore, treat as incorrect/correct) - • Comprehensive performance benchmarking suite - • Memory-efficient implementation with excellent scaling - • Production-ready with extensive test coverage + • Joint maximum likelihood by gradient ascent with an adaptive learning rate + • Convergence reporting, scale identification, and exclusion of extreme response patterns + • Missing data strategies (ignore, treat as incorrect/correct) + • Performance and convergence benchmark scripts - Perfect for researchers, data scientists, and developers working with#{" "} + Intended for researchers, data scientists, and developers working with#{" "} educational assessments, psychological measurements, or any binary response data where item and person parameters need to be estimated simultaneously. DESC diff --git a/lib/irt_ruby/estimation.rb b/lib/irt_ruby/estimation.rb new file mode 100644 index 0000000..659a6cf --- /dev/null +++ b/lib/irt_ruby/estimation.rb @@ -0,0 +1,95 @@ +# frozen_string_literal: true + +module IrtRuby + # Numerical helpers shared by IRT model implementations. + # @api private + module Estimation + # Learning-rate multiplier after an accepted step, so the rate can recover from earlier decays. + LEARNING_RATE_GROWTH = 1.1 + + module_function + + def log_sigmoid(x) + x >= 0 ? -Math.log(1.0 + Math.exp(-x)) : x - Math.log(1.0 + Math.exp(x)) + end + + # Returns [persons, items] flags marking which rows and columns have finite + # maximum-likelihood estimates. A person or item whose observed responses are + # all correct, all incorrect, or absent has none, so it is excluded; exclusion + # repeats until stable because dropping an item can make a person extreme. + def estimable(responses) + persons = Array.new(responses.size, true) + items = Array.new(responses.first.size, true) + + loop do + next_persons = persons.each_index.map do |i| + persons[i] && mixed?(items.each_index.select { |j| items[j] }.map { |j| responses[i][j] }) + end + next_items = items.each_index.map do |j| + items[j] && mixed?(next_persons.each_index.select { |i| next_persons[i] }.map { |i| responses[i][j] }) + end + return [persons, items] if next_persons == persons && next_items == items + + persons = next_persons + items = next_items + end + end + + def mixed?(values) + values.include?(0) && values.include?(1) + end + + def each_response(responses, persons, items) + responses.each_with_index do |row, i| + next unless persons[i] + + row.each_with_index do |value, j| + next if value.nil? || !items[j] + + yield i, j, value + end + end + end + + def response_counts(responses, persons, items) + person_counts = Array.new(responses.size, 0) + item_counts = Array.new(responses.first.size, 0) + each_response(responses, persons, items) do |i, j, _value| + person_counts[i] += 1 + item_counts[j] += 1 + end + [person_counts, item_counts] + end + + # True when every gradient component, averaged over the responses it sums, is below tolerance. + def stationary?(person_gradient, item_gradients, person_counts, item_counts, tolerance) + small_residuals?(person_gradient, person_counts, tolerance) && + item_gradients.all? { |gradient| small_residuals?(gradient, item_counts, tolerance) } + end + + def small_residuals?(gradient, counts, tolerance) + gradient.each_index.all? { |k| counts[k].zero? || (gradient[k] / counts[k]).abs < tolerance } + end + + # Drops gradient components that push a parameter further past the bound it sits on. + def project(gradient, values, min, max) + gradient.each_index.map do |k| + blocked = (values[k] <= min && gradient[k].negative?) || (values[k] >= max && gradient[k].positive?) + blocked ? 0.0 : gradient[k] + end + end + + def mean_and_sd(values) + mean = values.sum / values.size + [mean, Math.sqrt(values.sum { |v| (v - mean)**2 } / values.size)] + end + + def changes(current, previous, active) + current.each_index.select { |k| active[k] }.map { |k| (current[k] - previous[k]).abs } + end + + def estimates(values, active) + values.each_index.map { |k| active[k] ? values[k] : nil } + end + end +end diff --git a/lib/irt_ruby/model_options_validator.rb b/lib/irt_ruby/model_options_validator.rb index 56973b2..2936e49 100644 --- a/lib/irt_ruby/model_options_validator.rb +++ b/lib/irt_ruby/model_options_validator.rb @@ -5,10 +5,11 @@ module IrtRuby module ModelOptionsValidator module_function - def validate!(max_iter:, tolerance:, param_tolerance:, learning_rate:, decay_factor:) + def validate!(max_iter:, tolerance:, param_tolerance:, gradient_tolerance:, learning_rate:, decay_factor:) validate_positive_integer!(:max_iter, max_iter) validate_positive_finite_numeric!(:tolerance, tolerance) validate_positive_finite_numeric!(:param_tolerance, param_tolerance) + validate_positive_finite_numeric!(:gradient_tolerance, gradient_tolerance) validate_positive_finite_numeric!(:learning_rate, learning_rate) validate_decay_factor!(decay_factor) end diff --git a/lib/irt_ruby/rasch_model.rb b/lib/irt_ruby/rasch_model.rb index c6f44ff..851c2ec 100644 --- a/lib/irt_ruby/rasch_model.rb +++ b/lib/irt_ruby/rasch_model.rb @@ -1,5 +1,6 @@ # frozen_string_literal: true +require "irt_ruby/estimation" require "irt_ruby/model_options_validator" require "irt_ruby/response_data_validator" @@ -9,6 +10,8 @@ module IrtRuby # - Adaptive learning rate # - Missing data handling (skip nil) # - Multiple convergence checks (log-likelihood + parameter updates) + # - Exclusion of extreme response patterns, which have no finite estimates + # - Identification by centering abilities at zero class RaschModel MISSING_STRATEGIES = %i[ignore treat_as_incorrect treat_as_correct].freeze @@ -16,6 +19,7 @@ def initialize(data, max_iter: 1000, tolerance: 1e-6, param_tolerance: 1e-6, + gradient_tolerance: 1e-3, learning_rate: 0.01, decay_factor: 0.5, missing_strategy: :ignore, @@ -26,28 +30,33 @@ def initialize(data, ModelOptionsValidator.validate!(max_iter: max_iter, tolerance: tolerance, param_tolerance: param_tolerance, + gradient_tolerance: gradient_tolerance, learning_rate: learning_rate, decay_factor: decay_factor) - @data = data - @data_array = ResponseDataValidator.validate!(data) - num_rows = @data_array.size - num_cols = @data_array.first.size + data_array = ResponseDataValidator.validate!(data) + num_rows = data_array.size + num_cols = data_array.first.size raise ArgumentError, "missing_strategy must be one of #{MISSING_STRATEGIES}" unless MISSING_STRATEGIES.include?(missing_strategy) @missing_strategy = missing_strategy @random = seed.nil? ? nil : Random.new(seed) + @responses = data_array.map { |row| row.map { |resp| resolve_missing(resp).first } } + @active_persons, @active_items = Estimation.estimable(@responses) + @person_counts, @item_counts = Estimation.response_counts(@responses, @active_persons, @active_items) + # Initialize parameters near zero @abilities = Array.new(num_rows) { random_between(-0.25..0.25) } @difficulties = Array.new(num_cols) { random_between(-0.25..0.25) } - @max_iter = max_iter - @tolerance = tolerance - @param_tolerance = param_tolerance - @learning_rate = learning_rate - @decay_factor = decay_factor + @max_iter = max_iter + @tolerance = tolerance + @param_tolerance = param_tolerance + @gradient_tolerance = gradient_tolerance + @learning_rate = learning_rate + @decay_factor = decay_factor end def sigmoid(x) @@ -74,18 +83,9 @@ def resolve_missing(resp) def log_likelihood total_ll = 0.0 - @data_array.each_with_index do |row, i| - row.each_with_index do |resp, j| - value, skip = resolve_missing(resp) - next if skip - - prob = sigmoid(@abilities[i] - @difficulties[j]) - total_ll += if value == 1 - Math.log(prob + 1e-15) - else - Math.log((1 - prob) + 1e-15) - end - end + Estimation.each_response(@responses, @active_persons, @active_items) do |i, j, value| + z = @abilities[i] - @difficulties[j] + total_ll += Estimation.log_sigmoid(value == 1 ? z : -z) end total_ll end @@ -94,17 +94,11 @@ def compute_gradient grad_abilities = Array.new(@abilities.size, 0.0) grad_difficulties = Array.new(@difficulties.size, 0.0) - @data_array.each_with_index do |row, i| - row.each_with_index do |resp, j| - value, skip = resolve_missing(resp) - next if skip + Estimation.each_response(@responses, @active_persons, @active_items) do |i, j, value| + error = value - sigmoid(@abilities[i] - @difficulties[j]) - prob = sigmoid(@abilities[i] - @difficulties[j]) - error = value - prob - - grad_abilities[i] += error - grad_difficulties[j] -= error - end + grad_abilities[i] += error + grad_difficulties[j] -= error end [grad_abilities, grad_difficulties] @@ -126,40 +120,67 @@ def apply_gradient_update(grad_abilities, grad_difficulties) end def average_param_update(old_abilities, old_difficulties) - deltas = [] - @abilities.each_with_index do |a, i| - deltas << (a - old_abilities[i]).abs - end - @difficulties.each_with_index do |d, j| - deltas << (d - old_difficulties[j]).abs - end + deltas = Estimation.changes(@abilities, old_abilities, @active_persons) + + Estimation.changes(@difficulties, old_difficulties, @active_items) deltas.sum / deltas.size end + # Returns estimates with nil for persons and items excluded as extreme or + # unobserved, plus whether the convergence criteria were met. def fit - prev_ll = log_likelihood - - @max_iter.times do - grad_abilities, grad_difficulties = compute_gradient - - old_a, old_d = apply_gradient_update(grad_abilities, grad_difficulties) + converged = false + iterations = 0 + + if @active_persons.any? + prev_ll = log_likelihood + + @max_iter.times do + iterations += 1 + grad_abilities, grad_difficulties = compute_gradient + + old_a, old_d = apply_gradient_update(grad_abilities, grad_difficulties) + + current_ll = log_likelihood + param_delta = average_param_update(old_a, old_d) + + if current_ll < prev_ll + @abilities = old_a + @difficulties = old_d + @learning_rate *= @decay_factor + else + ll_diff = (current_ll - prev_ll).abs + if ll_diff < @tolerance && param_delta < @param_tolerance + converged = stationary? + break + end + + prev_ll = current_ll + @learning_rate *= Estimation::LEARNING_RATE_GROWTH + end + end - current_ll = log_likelihood - param_delta = average_param_update(old_a, old_d) + center_abilities + end - if current_ll < prev_ll - @abilities = old_a - @difficulties = old_d - @learning_rate *= @decay_factor - else - ll_diff = (current_ll - prev_ll).abs - break if ll_diff < @tolerance && param_delta < @param_tolerance + { + abilities: Estimation.estimates(@abilities, @active_persons), + difficulties: Estimation.estimates(@difficulties, @active_items), + converged: converged, + iterations: iterations + } + end - prev_ll = current_ll - end - end + def stationary? + grad_abilities, grad_difficulties = compute_gradient + Estimation.stationary?(grad_abilities, [grad_difficulties], @person_counts, @item_counts, @gradient_tolerance) + end + private :stationary? - { abilities: @abilities, difficulties: @difficulties } + def center_abilities + mean, = Estimation.mean_and_sd(@abilities.select.with_index { |_, i| @active_persons[i] }) + @abilities.map! { |theta| theta - mean } + @difficulties.map! { |b| b - mean } end + private :center_abilities end end diff --git a/lib/irt_ruby/three_parameter_model.rb b/lib/irt_ruby/three_parameter_model.rb index 28a58c3..7fccde1 100644 --- a/lib/irt_ruby/three_parameter_model.rb +++ b/lib/irt_ruby/three_parameter_model.rb @@ -1,5 +1,6 @@ # frozen_string_literal: true +require "irt_ruby/estimation" require "irt_ruby/model_options_validator" require "irt_ruby/response_data_validator" @@ -11,13 +12,18 @@ module IrtRuby # - Parameter clamping for discrimination, guessing # - Multiple convergence checks # - Separate gradient calculation & updates + # - Exclusion of extreme response patterns, which have no finite estimates + # - Identification by standardizing abilities to mean 0 and SD 1 class ThreeParameterModel MISSING_STRATEGIES = %i[ignore treat_as_incorrect treat_as_correct].freeze + DISCRIMINATION_BOUNDS = [0.01, 5.0].freeze + GUESSING_BOUNDS = [0.0, 0.35].freeze def initialize(data, max_iter: 1000, tolerance: 1e-6, param_tolerance: 1e-6, + gradient_tolerance: 1e-3, learning_rate: 0.01, decay_factor: 0.5, missing_strategy: :ignore, @@ -25,30 +31,35 @@ def initialize(data, ModelOptionsValidator.validate!(max_iter: max_iter, tolerance: tolerance, param_tolerance: param_tolerance, + gradient_tolerance: gradient_tolerance, learning_rate: learning_rate, decay_factor: decay_factor) - @data = data - @data_array = ResponseDataValidator.validate!(data) - num_rows = @data_array.size - num_cols = @data_array.first.size + data_array = ResponseDataValidator.validate!(data) + num_rows = data_array.size + num_cols = data_array.first.size raise ArgumentError, "missing_strategy must be one of #{MISSING_STRATEGIES}" unless MISSING_STRATEGIES.include?(missing_strategy) @missing_strategy = missing_strategy @random = seed.nil? ? nil : Random.new(seed) + @responses = data_array.map { |row| row.map { |resp| resolve_missing(resp).first } } + @active_persons, @active_items = Estimation.estimable(@responses) + @person_counts, @item_counts = Estimation.response_counts(@responses, @active_persons, @active_items) + # Initialize parameters @abilities = Array.new(num_rows) { random_between(-0.25..0.25) } @difficulties = Array.new(num_cols) { random_between(-0.25..0.25) } @discriminations = Array.new(num_cols) { random_between(0.5..1.5) } @guessings = Array.new(num_cols) { random_between(0.0..0.3) } - @max_iter = max_iter - @tolerance = tolerance - @param_tolerance = param_tolerance - @learning_rate = learning_rate - @decay_factor = decay_factor + @max_iter = max_iter + @tolerance = tolerance + @param_tolerance = param_tolerance + @gradient_tolerance = gradient_tolerance + @learning_rate = learning_rate + @decay_factor = decay_factor end def sigmoid(x) @@ -80,22 +91,15 @@ def resolve_missing(resp) def log_likelihood ll = 0.0 - @data_array.each_with_index do |row, i| - row.each_with_index do |resp, j| - value, skip = resolve_missing(resp) - next if skip - - prob = probability(@abilities[i], - @discriminations[j], - @difficulties[j], - @guessings[j]) - - ll += if value == 1 - Math.log(prob + 1e-15) - else - Math.log((1 - prob) + 1e-15) - end - end + Estimation.each_response(@responses, @active_persons, @active_items) do |i, j, value| + z = @discriminations[j] * (@abilities[i] - @difficulties[j]) + c = @guessings[j] + + ll += if value == 1 + c.zero? ? Estimation.log_sigmoid(z) : Math.log(c + ((1.0 - c) * sigmoid(z))) + else + Math.log(1.0 - c) + Estimation.log_sigmoid(-z) + end end ll end @@ -106,25 +110,23 @@ def compute_gradient grad_discriminations = Array.new(@discriminations.size, 0.0) grad_guessings = Array.new(@guessings.size, 0.0) - @data_array.each_with_index do |row, i| - row.each_with_index do |resp, j| - value, skip = resolve_missing(resp) - next if skip - - theta = @abilities[i] - a = @discriminations[j] - b = @difficulties[j] - c = @guessings[j] - - prob = probability(theta, a, b, c) - error = value - prob - - grad_abilities[i] += error * a * (1 - c) - grad_difficulties[j] -= error * a * (1 - c) - grad_discriminations[j] += error * (theta - b) * (1 - c) - - grad_guessings[j] += error * 1.0 - end + Estimation.each_response(@responses, @active_persons, @active_items) do |i, j, value| + theta = @abilities[i] + a = @discriminations[j] + b = @difficulties[j] + c = @guessings[j] + + # With P = c + (1 - c) * s and s = sigmoid(a * (theta - b)): + # dLL/dz = (u - P) * s / P and dLL/dc = (u - P) / ((1 - c) * P). + s = sigmoid(a * (theta - b)) + prob = c + ((1.0 - c) * s) + error = value - prob + dz = prob.zero? ? error : error * s / prob + + grad_abilities[i] += dz * a + grad_difficulties[j] -= dz * a + grad_discriminations[j] += dz * (theta - b) + grad_guessings[j] += prob.zero? ? 0.0 : error / ((1.0 - c) * prob) end [grad_abilities, grad_difficulties, grad_discriminations, grad_guessings] @@ -145,59 +147,89 @@ def apply_gradient_update(ga, gd, gdisc, gc) end @discriminations.each_index do |j| - @discriminations[j] += @learning_rate * gdisc[j] - @discriminations[j] = 0.01 if @discriminations[j] < 0.01 - @discriminations[j] = 5.0 if @discriminations[j] > 5.0 + @discriminations[j] = (@discriminations[j] + (@learning_rate * gdisc[j])).clamp(*DISCRIMINATION_BOUNDS) end @guessings.each_index do |j| - @guessings[j] += @learning_rate * gc[j] - @guessings[j] = 0.0 if @guessings[j] < 0.0 - @guessings[j] = 0.35 if @guessings[j] > 0.35 + @guessings[j] = (@guessings[j] + (@learning_rate * gc[j])).clamp(*GUESSING_BOUNDS) end [old_a, old_d, old_disc, old_c] end def average_param_update(old_a, old_d, old_disc, old_c) - deltas = [] - @abilities.each_with_index { |x, i| deltas << (x - old_a[i]).abs } - @difficulties.each_with_index { |x, j| deltas << (x - old_d[j]).abs } - @discriminations.each_with_index { |x, j| deltas << (x - old_disc[j]).abs } - @guessings.each_with_index { |x, j| deltas << (x - old_c[j]).abs } + deltas = Estimation.changes(@abilities, old_a, @active_persons) + + Estimation.changes(@difficulties, old_d, @active_items) + + Estimation.changes(@discriminations, old_disc, @active_items) + + Estimation.changes(@guessings, old_c, @active_items) deltas.sum / deltas.size end + # Returns estimates with nil for persons and items excluded as extreme or + # unobserved, plus whether the convergence criteria were met. def fit - prev_ll = log_likelihood - - @max_iter.times do - ga, gd, gdisc, gc = compute_gradient - old_a, old_d, old_disc, old_c = apply_gradient_update(ga, gd, gdisc, gc) - - curr_ll = log_likelihood - param_delta = average_param_update(old_a, old_d, old_disc, old_c) - - if curr_ll < prev_ll - @abilities = old_a - @difficulties = old_d - @discriminations = old_disc - @guessings = old_c - @learning_rate *= @decay_factor - else - ll_diff = (curr_ll - prev_ll).abs - break if ll_diff < @tolerance && param_delta < @param_tolerance - - prev_ll = curr_ll + converged = false + iterations = 0 + + if @active_persons.any? + prev_ll = log_likelihood + + @max_iter.times do + iterations += 1 + ga, gd, gdisc, gc = compute_gradient + old_a, old_d, old_disc, old_c = apply_gradient_update(ga, gd, gdisc, gc) + + curr_ll = log_likelihood + param_delta = average_param_update(old_a, old_d, old_disc, old_c) + + if curr_ll < prev_ll + @abilities = old_a + @difficulties = old_d + @discriminations = old_disc + @guessings = old_c + @learning_rate *= @decay_factor + else + ll_diff = (curr_ll - prev_ll).abs + if ll_diff < @tolerance && param_delta < @param_tolerance + converged = stationary? + break + end + + prev_ll = curr_ll + @learning_rate *= Estimation::LEARNING_RATE_GROWTH + end end + + standardize_abilities end { - abilities: @abilities, - difficulties: @difficulties, - discriminations: @discriminations, - guessings: @guessings + abilities: Estimation.estimates(@abilities, @active_persons), + difficulties: Estimation.estimates(@difficulties, @active_items), + discriminations: Estimation.estimates(@discriminations, @active_items), + guessings: Estimation.estimates(@guessings, @active_items), + converged: converged, + iterations: iterations } end + + def stationary? + ga, gd, gdisc, gc = compute_gradient + gdisc = Estimation.project(gdisc, @discriminations, *DISCRIMINATION_BOUNDS) + gc = Estimation.project(gc, @guessings, *GUESSING_BOUNDS) + Estimation.stationary?(ga, [gd, gdisc, gc], @person_counts, @item_counts, @gradient_tolerance) + end + private :stationary? + + # Rescales the latent trait without changing the likelihood: theta' = (theta - m) / s, + # b' = (b - m) / s, a' = a * s. + def standardize_abilities + mean, sd = Estimation.mean_and_sd(@abilities.select.with_index { |_, i| @active_persons[i] }) + sd = 1.0 if sd.zero? + @abilities.map! { |theta| (theta - mean) / sd } + @difficulties.map! { |b| (b - mean) / sd } + @discriminations.map! { |a| a * sd } + end + private :standardize_abilities end end diff --git a/lib/irt_ruby/two_parameter_model.rb b/lib/irt_ruby/two_parameter_model.rb index 1491fb5..a4058a4 100644 --- a/lib/irt_ruby/two_parameter_model.rb +++ b/lib/irt_ruby/two_parameter_model.rb @@ -1,5 +1,6 @@ # frozen_string_literal: true +require "irt_ruby/estimation" require "irt_ruby/model_options_validator" require "irt_ruby/response_data_validator" @@ -11,39 +12,48 @@ module IrtRuby # - Parameter clamping for discrimination # - Multiple convergence checks # - Separate gradient calculation & parameter update + # - Exclusion of extreme response patterns, which have no finite estimates + # - Identification by standardizing abilities to mean 0 and SD 1 class TwoParameterModel MISSING_STRATEGIES = %i[ignore treat_as_incorrect treat_as_correct].freeze + DISCRIMINATION_BOUNDS = [0.01, 5.0].freeze def initialize(data, max_iter: 1000, tolerance: 1e-6, param_tolerance: 1e-6, + gradient_tolerance: 1e-3, learning_rate: 0.01, decay_factor: 0.5, missing_strategy: :ignore, seed: nil) ModelOptionsValidator.validate!(max_iter: max_iter, tolerance: tolerance, param_tolerance: param_tolerance, + gradient_tolerance: gradient_tolerance, learning_rate: learning_rate, decay_factor: decay_factor) - @data = data - @data_array = ResponseDataValidator.validate!(data) - num_rows = @data_array.size - num_cols = @data_array.first.size + data_array = ResponseDataValidator.validate!(data) + num_rows = data_array.size + num_cols = data_array.first.size raise ArgumentError, "missing_strategy must be one of #{MISSING_STRATEGIES}" unless MISSING_STRATEGIES.include?(missing_strategy) @missing_strategy = missing_strategy @random = seed.nil? ? nil : Random.new(seed) + @responses = data_array.map { |row| row.map { |resp| resolve_missing(resp).first } } + @active_persons, @active_items = Estimation.estimable(@responses) + @person_counts, @item_counts = Estimation.response_counts(@responses, @active_persons, @active_items) + # Initialize parameters # Typically: ability ~ 0, difficulty ~ 0, discrimination ~ 1 @abilities = Array.new(num_rows) { random_between(-0.25..0.25) } @difficulties = Array.new(num_cols) { random_between(-0.25..0.25) } @discriminations = Array.new(num_cols) { random_between(0.5..1.5) } - @max_iter = max_iter - @tolerance = tolerance - @param_tolerance = param_tolerance - @learning_rate = learning_rate - @decay_factor = decay_factor + @max_iter = max_iter + @tolerance = tolerance + @param_tolerance = param_tolerance + @gradient_tolerance = gradient_tolerance + @learning_rate = learning_rate + @decay_factor = decay_factor end def sigmoid(x) @@ -70,18 +80,9 @@ def resolve_missing(resp) def log_likelihood ll = 0.0 - @data_array.each_with_index do |row, i| - row.each_with_index do |resp, j| - value, skip = resolve_missing(resp) - next if skip - - prob = sigmoid(@discriminations[j] * (@abilities[i] - @difficulties[j])) - ll += if value == 1 - Math.log(prob + 1e-15) - else - Math.log((1 - prob) + 1e-15) - end - end + Estimation.each_response(@responses, @active_persons, @active_items) do |i, j, value| + z = @discriminations[j] * (@abilities[i] - @difficulties[j]) + ll += Estimation.log_sigmoid(value == 1 ? z : -z) end ll end @@ -91,18 +92,13 @@ def compute_gradient grad_difficulties = Array.new(@difficulties.size, 0.0) grad_discriminations = Array.new(@discriminations.size, 0.0) - @data_array.each_with_index do |row, i| - row.each_with_index do |resp, j| - value, skip = resolve_missing(resp) - next if skip - - prob = sigmoid(@discriminations[j] * (@abilities[i] - @difficulties[j])) - error = value - prob + Estimation.each_response(@responses, @active_persons, @active_items) do |i, j, value| + prob = sigmoid(@discriminations[j] * (@abilities[i] - @difficulties[j])) + error = value - prob - grad_abilities[i] += error * @discriminations[j] - grad_difficulties[j] -= error * @discriminations[j] - grad_discriminations[j] += error * (@abilities[i] - @difficulties[j]) - end + grad_abilities[i] += error * @discriminations[j] + grad_difficulties[j] -= error * @discriminations[j] + grad_discriminations[j] += error * (@abilities[i] - @difficulties[j]) end [grad_abilities, grad_difficulties, grad_discriminations] @@ -122,50 +118,81 @@ def apply_gradient_update(ga, gd, gdisc) end @discriminations.each_index do |j| - @discriminations[j] += @learning_rate * gdisc[j] - @discriminations[j] = 0.01 if @discriminations[j] < 0.01 - @discriminations[j] = 5.0 if @discriminations[j] > 5.0 + @discriminations[j] = (@discriminations[j] + (@learning_rate * gdisc[j])).clamp(*DISCRIMINATION_BOUNDS) end [old_a, old_d, old_disc] end def average_param_update(old_a, old_d, old_disc) - deltas = [] - @abilities.each_with_index { |x, i| deltas << (x - old_a[i]).abs } - @difficulties.each_with_index { |x, j| deltas << (x - old_d[j]).abs } - @discriminations.each_with_index { |x, j| deltas << (x - old_disc[j]).abs } + deltas = Estimation.changes(@abilities, old_a, @active_persons) + + Estimation.changes(@difficulties, old_d, @active_items) + + Estimation.changes(@discriminations, old_disc, @active_items) deltas.sum / deltas.size end + # Returns estimates with nil for persons and items excluded as extreme or + # unobserved, plus whether the convergence criteria were met. def fit - prev_ll = log_likelihood - - @max_iter.times do - ga, gd, gdisc = compute_gradient - old_a, old_d, old_disc = apply_gradient_update(ga, gd, gdisc) - - curr_ll = log_likelihood - param_delta = average_param_update(old_a, old_d, old_disc) - - if curr_ll < prev_ll - @abilities = old_a - @difficulties = old_d - @discriminations = old_disc - @learning_rate *= @decay_factor - else - ll_diff = (curr_ll - prev_ll).abs - break if ll_diff < @tolerance && param_delta < @param_tolerance - - prev_ll = curr_ll + converged = false + iterations = 0 + + if @active_persons.any? + prev_ll = log_likelihood + + @max_iter.times do + iterations += 1 + ga, gd, gdisc = compute_gradient + old_a, old_d, old_disc = apply_gradient_update(ga, gd, gdisc) + + curr_ll = log_likelihood + param_delta = average_param_update(old_a, old_d, old_disc) + + if curr_ll < prev_ll + @abilities = old_a + @difficulties = old_d + @discriminations = old_disc + @learning_rate *= @decay_factor + else + ll_diff = (curr_ll - prev_ll).abs + if ll_diff < @tolerance && param_delta < @param_tolerance + converged = stationary? + break + end + + prev_ll = curr_ll + @learning_rate *= Estimation::LEARNING_RATE_GROWTH + end end + + standardize_abilities end { - abilities: @abilities, - difficulties: @difficulties, - discriminations: @discriminations + abilities: Estimation.estimates(@abilities, @active_persons), + difficulties: Estimation.estimates(@difficulties, @active_items), + discriminations: Estimation.estimates(@discriminations, @active_items), + converged: converged, + iterations: iterations } end + + def stationary? + ga, gd, gdisc = compute_gradient + gdisc = Estimation.project(gdisc, @discriminations, *DISCRIMINATION_BOUNDS) + Estimation.stationary?(ga, [gd, gdisc], @person_counts, @item_counts, @gradient_tolerance) + end + private :stationary? + + # Rescales the latent trait without changing the likelihood: theta' = (theta - m) / s, + # b' = (b - m) / s, a' = a * s. + def standardize_abilities + mean, sd = Estimation.mean_and_sd(@abilities.select.with_index { |_, i| @active_persons[i] }) + sd = 1.0 if sd.zero? + @abilities.map! { |theta| (theta - mean) / sd } + @difficulties.map! { |b| (b - mean) / sd } + @discriminations.map! { |a| a * sd } + end + private :standardize_abilities end end diff --git a/spec/irt_ruby/estimation_spec.rb b/spec/irt_ruby/estimation_spec.rb new file mode 100644 index 0000000..3d1e915 --- /dev/null +++ b/spec/irt_ruby/estimation_spec.rb @@ -0,0 +1,174 @@ +# frozen_string_literal: true + +require "spec_helper" + +RSpec.describe "Estimation correctness" do + models = { + IrtRuby::RaschModel => %i[abilities difficulties], + IrtRuby::TwoParameterModel => %i[abilities difficulties discriminations], + IrtRuby::ThreeParameterModel => %i[abilities difficulties discriminations guessings] + } + + def numeric_derivative(model, parameter, index, step = 1e-6) + values = model.instance_variable_get("@#{parameter}") + original = values[index] + values[index] = original + step + upper = model.log_likelihood + values[index] = original - step + lower = model.log_likelihood + values[index] = original + (upper - lower) / (2 * step) + end + + def correlation(xs, ys) + mean_x = xs.sum / xs.size + mean_y = ys.sum / ys.size + covariance = xs.zip(ys).sum { |x, y| (x - mean_x) * (y - mean_y) } + covariance / Math.sqrt(xs.sum { |x| (x - mean_x)**2 } * ys.sum { |y| (y - mean_y)**2 }) + end + + def simulate(rng, abilities, difficulties, discriminations, guessings) + abilities.map do |theta| + difficulties.each_index.map do |j| + prob = guessings[j] + ((1 - guessings[j]) / (1 + Math.exp(-discriminations[j] * (theta - difficulties[j])))) + rng.rand < prob ? 1 : 0 + end + end + end + + let(:mixed_data) do + rng = Random.new(7) + Array.new(12) { Array.new(6) { [0, 1, 1, nil][rng.rand(4)] } } + end + + models.each do |model_class, parameters| + describe model_class do + %i[ignore treat_as_incorrect treat_as_correct].each do |strategy| + it "has analytic gradients matching finite differences with #{strategy}" do + model = model_class.new(mixed_data, missing_strategy: strategy, seed: 3) + gradients = model.compute_gradient + + parameters.each_with_index do |parameter, n| + gradients[n].each_index do |k| + next if gradients[n][k].zero? && numeric_derivative(model, parameter, k).zero? + + expect(gradients[n][k]).to be_within(1e-5).of(numeric_derivative(model, parameter, k)) + end + end + end + end + + it "returns nil for persons and items with extreme or unobserved responses" do + data = [ + [1, 0, 1, 1, nil], + [0, 1, 0, 1, nil], + [1, 1, 0, 1, nil], + [1, 1, 1, 1, nil], + [0, 0, 0, 1, nil], + [nil, nil, nil, nil, nil] + ] + results = model_class.new(data, seed: 1).fit + + expect(results[:abilities].map(&:nil?)).to eq([false, false, false, true, true, true]) + expect(results[:difficulties].map(&:nil?)).to eq([false, false, false, true, true]) + end + + it "repeats exclusion when dropping items makes a person extreme" do + data = [ + [1, 1, 0], + [1, 0, 0], + [0, 1, 0] + ] + results = model_class.new(data, seed: 1).fit + + expect(results[:difficulties].map(&:nil?)).to eq([false, false, true]) + expect(results[:abilities].map(&:nil?)).to eq([true, false, false]) + end + + it "fixes the latent scale of the returned abilities" do + rng = Random.new(11) + data = simulate(rng, Array.new(60) { (rng.rand * 4) - 2 }, Array.new(8) { (rng.rand * 2) - 1 }, + Array.new(8, 1.0), Array.new(8, 0.0)) + abilities = model_class.new(data, max_iter: 50, seed: 1).fit[:abilities].compact + mean, sd = IrtRuby::Estimation.mean_and_sd(abilities) + + expect(mean).to be_within(1e-12).of(0.0) + expect(sd).to be_within(1e-12).of(1.0) unless model_class == IrtRuby::RaschModel + end + + it "keeps the likelihood unchanged when fixing the scale" do + model = model_class.new(mixed_data, max_iter: 1, seed: 1) + model.instance_variable_set(:@max_iter, 0) + before = model.log_likelihood + model.fit + + expect(model.log_likelihood).to be_within(1e-9).of(before) + end + + it "does not report convergence when the learning rate is too small to move" do + results = model_class.new(mixed_data, learning_rate: 1e-20, seed: 1).fit + + expect(results[:converged]).to be(false) + expect(results[:iterations]).to eq(1) + end + + it "reports exhausting max_iter as not converged" do + results = model_class.new(mixed_data, max_iter: 3, seed: 1).fit + + expect(results).to include(converged: false, iterations: 3) + end + + it "evaluates badly mispredicted responses without saturating the likelihood" do + model = model_class.new([[0, 1], [1, 0]], seed: 1) + model.instance_variable_set(:@abilities, [40.0, 0.0]) + model.instance_variable_set(:@difficulties, [0.0, 0.0]) + model.instance_variable_set(:@discriminations, [1.0, 1.0]) if parameters.include?(:discriminations) + model.instance_variable_set(:@guessings, [0.0, 0.0]) if parameters.include?(:guessings) + + expect(model.log_likelihood).to be_within(1e-9).of(-40.0 - (2 * Math.log(2))) + end + + it "is not affected by later changes to the input or the returned arrays" do + data = [[1, 0], [0, 1]] + model = model_class.new(data, max_iter: 5, seed: 1) + data[0] << 1 + results = model.fit + before = model.log_likelihood + results[:abilities][0] = 100.0 + + expect(model.log_likelihood).to eq(before) + end + end + end + + it "reaches the exact Rasch joint estimates on the README data" do + results = IrtRuby::RaschModel.new([[1, 0, 1], [0, 1, 0], [1, 1, 1]], learning_rate: 0.5, seed: 1).fit + + expect(results[:converged]).to be(true) + expect(results[:abilities][0]).to be_within(1e-3).of(Math.log(2)) + expect(results[:abilities][1]).to be_within(1e-3).of(-Math.log(2)) + expect(results[:abilities][2]).to be_nil + results[:difficulties].each { |b| expect(b).to be_within(1e-3).of(0.0) } + end + + it "recovers simulated Rasch difficulties" do + rng = Random.new(42) + difficulties = Array.new(10) { (rng.rand * 3) - 1.5 } + data = simulate(rng, Array.new(300) { (rng.rand * 4) - 2 }, difficulties, Array.new(10, 1.0), Array.new(10, 0.0)) + results = IrtRuby::RaschModel.new(data, learning_rate: 0.05, seed: 1).fit + + expect(results[:converged]).to be(true) + expect(correlation(difficulties, results[:difficulties])).to be > 0.97 + end + + it "recovers simulated 2PL difficulties and discriminations" do + rng = Random.new(42) + difficulties = Array.new(10) { (rng.rand * 3) - 1.5 } + discriminations = Array.new(10) { 0.6 + (rng.rand * 1.4) } + data = simulate(rng, Array.new(300) { (rng.rand * 4) - 2 }, difficulties, discriminations, Array.new(10, 0.0)) + results = IrtRuby::TwoParameterModel.new(data, learning_rate: 0.05, max_iter: 500, seed: 1).fit + + expect(correlation(difficulties, results[:difficulties])).to be > 0.95 + expect(correlation(discriminations, results[:discriminations])).to be > 0.7 + end +end diff --git a/spec/irt_ruby/three_parameter_model_spec.rb b/spec/irt_ruby/three_parameter_model_spec.rb index a0b6178..af704a2 100644 --- a/spec/irt_ruby/three_parameter_model_spec.rb +++ b/spec/irt_ruby/three_parameter_model_spec.rb @@ -245,22 +245,16 @@ end end - context "Known parameter test (optional)" do - it "attempts to recover small synthetic data parameters" do - data = [ - [1, 1], - [1, 1] - ] - - model = described_class.new(data, max_iter: 200, learning_rate: 0.05) + context "Extreme response patterns" do + it "returns nil for every parameter when all responses are correct" do + model = described_class.new([[1, 1], [1, 1]], max_iter: 200, learning_rate: 0.05) results = model.fit - results[:discriminations].each do |disc| - expect(disc).to be_between(0.01, 5.0) - end - results[:guessings].each do |g| - expect(g).to be_between(0.0, 0.35) - end + expect(results[:abilities]).to eq([nil, nil]) + expect(results[:difficulties]).to eq([nil, nil]) + expect(results[:discriminations]).to eq([nil, nil]) + expect(results[:guessings]).to eq([nil, nil]) + expect(results[:iterations]).to eq(0) end end end diff --git a/spec/irt_ruby/two_parameter_model_spec.rb b/spec/irt_ruby/two_parameter_model_spec.rb index fbc2a16..5bebece 100644 --- a/spec/irt_ruby/two_parameter_model_spec.rb +++ b/spec/irt_ruby/two_parameter_model_spec.rb @@ -202,18 +202,15 @@ end end - context "Known parameter test (optional)" do - it "checks parameter ranges on a small synthetic dataset" do - data = [ - [1, 1], - [1, 1] - ] - model = described_class.new(data, max_iter: 200, learning_rate: 0.05) + context "Extreme response patterns" do + it "returns nil for every parameter when all responses are correct" do + model = described_class.new([[1, 1], [1, 1]], max_iter: 200, learning_rate: 0.05) results = model.fit - results[:discriminations].each do |disc| - expect(disc).to be_between(0.01, 5.0) - end + expect(results[:abilities]).to eq([nil, nil]) + expect(results[:difficulties]).to eq([nil, nil]) + expect(results[:discriminations]).to eq([nil, nil]) + expect(results[:iterations]).to eq(0) end end end diff --git a/spec/spec_helper.rb b/spec/spec_helper.rb index 444abab..8e088f0 100644 --- a/spec/spec_helper.rb +++ b/spec/spec_helper.rb @@ -48,6 +48,7 @@ max_iter: [0, -1, 1.5, "100", nil], tolerance: [0, -1e-6, Float::INFINITY, -Float::INFINITY, Float::NAN, Complex(1, 0), "1e-6", nil], param_tolerance: [0, -1e-6, Float::INFINITY, -Float::INFINITY, Float::NAN, Complex(1, 0), "1e-6", nil], + gradient_tolerance: [0, -1e-6, Float::INFINITY, -Float::INFINITY, Float::NAN, Complex(1, 0), "1e-3", nil], learning_rate: [0, -0.01, Float::INFINITY, -Float::INFINITY, Float::NAN, Complex(0.01, 0), "0.01", nil], decay_factor: [0, 1, -0.1, 1.1, Float::INFINITY, -Float::INFINITY, Float::NAN, Complex(0.5, 0), "0.5", nil] }.each do |option, invalid_values| From 579f6cd4defee1d09a2f0231779ceb5d209210ca Mon Sep 17 00:00:00 2001 From: Alex Kholodniak Date: Sun, 11 Oct 2026 07:08:50 +0300 Subject: [PATCH 2/6] fix: estimate only the identified part of the response graph Mixed response rows and columns do not guarantee finite joint estimates: perfectly separated or disconnected groups still diverged or were left unidentified. Keep the largest strongly connected component of the response graph instead, which also covers extreme patterns. Standardizing the internal parameters pushed discriminations past their bounds and stalled later fit calls, and stationarity was judged on an arbitrary internal scale. Report rescaled copies and check stationarity on the reported scale. Under 3PL, all-incorrect respondents still inform guessing through log(1 - c), and the guessing gradient no longer drops to zero when the response function underflows. --- CHANGELOG.md | 14 ++--- README.md | 7 ++- lib/irt_ruby.rb | 2 + lib/irt_ruby/estimation.rb | 35 +++-------- lib/irt_ruby/rasch_model.rb | 22 +++---- lib/irt_ruby/response_graph.rb | 85 ++++++++++++++++++++++++++ lib/irt_ruby/three_parameter_model.rb | 67 ++++++++++++-------- lib/irt_ruby/two_parameter_model.rb | 35 +++++------ spec/irt_ruby/estimation_spec.rb | 88 +++++++++++++++++++++++++-- 9 files changed, 254 insertions(+), 101 deletions(-) create mode 100644 lib/irt_ruby/response_graph.rb diff --git a/CHANGELOG.md b/CHANGELOG.md index e5706cc..99985d9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,23 +6,23 @@ All notable changes to this project are documented in this file. ### Fixed - Fixed 3PL gradients for ability, difficulty, discrimination, and guessing; they now match finite differences of the log-likelihood. -- Respondents and items with all-correct, all-incorrect, or unobserved response patterns are excluded from estimation and returned as `nil` instead of drifting with `max_iter` or keeping random initial values. +- Only the largest strongly connected part of the response graph is estimated; other respondents and items are returned as `nil` instead of drifting with `max_iter` or keeping random initial values. This covers all-correct, all-incorrect, and unobserved rows or columns as well as perfectly separated or disconnected groups. +- 3PL keeps the limiting `log(1 - c)` contribution of excluded respondents who answered every estimated item incorrectly, and computes guessing gradients stably when the response function underflows. - Log-likelihoods are computed with a numerically stable log-sigmoid instead of adding `1e-15` inside `Math.log`, which flattened the objective for badly mispredicted responses. - Models copy response data on construction and return fresh arrays from `fit`, so later mutation of either does not change model state. ### Changed -- `fit` fixes the latent scale: Rasch abilities are centered at 0; 2PL and 3PL abilities are standardized to mean 0 and SD 1 with likelihood-preserving item transformations. -- `fit` results include `converged:` and `iterations:`. `converged` is true only at a stationary point, controlled by the new `gradient_tolerance:` option (default `1e-3`). +- `fit` reports results on a fixed latent scale: Rasch abilities are centered at 0; 2PL and 3PL abilities are standardized to mean 0 and SD 1 with likelihood-preserving item transformations. +- `fit` results include `converged:` and `iterations:`. `converged` is true only at a stationary point on the reported scale, controlled by the new `gradient_tolerance:` option (default `1e-3`). +- The learning rate grows by 10% after each accepted step, so it can recover from earlier decreases. - Clarified the gem summary and description to describe implemented capabilities rather than production-readiness claims. +- Added CI RuboCop and gem package build/install verification, and documented the local quality and packaging commands for contributors. +- Clarified that response data must be a `Matrix` or array of arrays containing only integer `0`, integer `1`, or `nil`; floats, strings, booleans, and other values are rejected. ### Added - Added finite-difference gradient, identification, extreme-pattern, and parameter-recovery specs. - Added executable specs for README usage examples to keep the documented Rasch, 2PL, 3PL, and missing-data flows aligned with the public API. -### Changed -- Added CI RuboCop and gem package build/install verification, and documented the local quality and packaging commands for contributors. -- Clarified that response data must be a `Matrix` or array of arrays containing only integer `0`, integer `1`, or `nil`; floats, strings, booleans, and other values are rejected. - --- ## [0.3.0] - 2025-06-13 diff --git a/README.md b/README.md index 33487d5..1fe2f14 100644 --- a/README.md +++ b/README.md @@ -155,14 +155,15 @@ For 2PL and 3PL: - **Discriminations** (`a`) are clamped between `0.01` and `5.0` during optimization. - **Guessings** (`c`, 3PL only) are clamped to `[0.0, 0.35]`. -These bounds are model restrictions: a 3PL item with a true lower asymptote above `0.35` cannot be represented. Because the final scale is standardized after optimization (see below), returned discriminations can fall slightly outside the optimization bounds. +These bounds are model restrictions: a 3PL item with a true lower asymptote above `0.35` cannot be represented. The bounds apply to the internal optimization scale; returned discriminations are multiplied by the internal ability SD (see below), so they can lie outside `[0.01, 5.0]`. ## Estimation Details All models estimate person and item parameters jointly by maximizing the likelihood (joint maximum likelihood, JMLE). -- **Extreme and unobserved patterns.** A respondent who answered every observed item correctly, or every item incorrectly, has no finite maximum-likelihood ability; the same holds for items that everyone answered correctly or incorrectly, and for rows or columns with no observed responses. These persons and items are excluded from estimation and returned as `nil`. Exclusion is repeated until stable, since removing an item can make a respondent's remaining pattern extreme. Missing-data strategies are applied first, so `:treat_as_correct` and `:treat_as_incorrect` can create extreme patterns. -- **Scale identification.** The latent scale is fixed after fitting: Rasch abilities are centered to mean 0; 2PL and 3PL abilities are standardized to mean 0 and standard deviation 1, with difficulties and discriminations transformed so that response probabilities are unchanged. Only estimable respondents enter these statistics. +- **Which parameters are estimated.** Think of the responses as a directed graph: a correct answer links the item to the respondent, an incorrect answer links the respondent to the item. Finite estimates on a common scale exist only within a strongly connected part of this graph, so IrtRuby estimates the largest such part and returns `nil` for every other respondent and item. This excludes respondents who answered every observed item correctly or every item incorrectly, items that everyone answered correctly or incorrectly, rows and columns with no observed responses, groups that are perfectly separated from the rest (for example, a group of respondents who answered every item of another group correctly), and groups that share no responses with the rest. Missing-data strategies are applied first, so `:treat_as_correct` and `:treat_as_incorrect` can create such patterns. +- **3PL and all-incorrect respondents.** Under 3PL, an excluded respondent who answered every estimated item incorrectly still informs the guessing parameters: as their ability tends to minus infinity each response has probability `1 - c`, and that limiting term stays in the likelihood. +- **Scale identification.** Results are reported on a fixed scale: Rasch abilities are centered to mean 0; 2PL and 3PL abilities are standardized to mean 0 and standard deviation 1, with difficulties and discriminations transformed so that response probabilities are unchanged. Only estimated respondents enter these statistics. The model keeps its internal parameters, so calling `fit` again continues from where it stopped; convergence is checked on the reported scale. - **Known limitations.** JMLE item estimates are biased when tests are short, and the bias does not vanish as the number of respondents grows. 2PL and 3PL JMLE is also slow to converge with plain gradient ascent. For high-stakes calibration, compare against marginal maximum likelihood software such as R's `mirt`. ## Performance Benchmarks diff --git a/lib/irt_ruby.rb b/lib/irt_ruby.rb index 6479b61..88ff502 100644 --- a/lib/irt_ruby.rb +++ b/lib/irt_ruby.rb @@ -4,6 +4,8 @@ require "matrix" require "irt_ruby/model_options_validator" require "irt_ruby/response_data_validator" +require "irt_ruby/estimation" +require "irt_ruby/response_graph" require "irt_ruby/rasch_model" require "irt_ruby/two_parameter_model" require "irt_ruby/three_parameter_model" diff --git a/lib/irt_ruby/estimation.rb b/lib/irt_ruby/estimation.rb index 659a6cf..1068a24 100644 --- a/lib/irt_ruby/estimation.rb +++ b/lib/irt_ruby/estimation.rb @@ -13,32 +13,6 @@ def log_sigmoid(x) x >= 0 ? -Math.log(1.0 + Math.exp(-x)) : x - Math.log(1.0 + Math.exp(x)) end - # Returns [persons, items] flags marking which rows and columns have finite - # maximum-likelihood estimates. A person or item whose observed responses are - # all correct, all incorrect, or absent has none, so it is excluded; exclusion - # repeats until stable because dropping an item can make a person extreme. - def estimable(responses) - persons = Array.new(responses.size, true) - items = Array.new(responses.first.size, true) - - loop do - next_persons = persons.each_index.map do |i| - persons[i] && mixed?(items.each_index.select { |j| items[j] }.map { |j| responses[i][j] }) - end - next_items = items.each_index.map do |j| - items[j] && mixed?(next_persons.each_index.select { |i| next_persons[i] }.map { |i| responses[i][j] }) - end - return [persons, items] if next_persons == persons && next_items == items - - persons = next_persons - items = next_items - end - end - - def mixed?(values) - values.include?(0) && values.include?(1) - end - def each_response(responses, persons, items) responses.each_with_index do |row, i| next unless persons[i] @@ -79,6 +53,15 @@ def project(gradient, values, min, max) end end + # Location and scale of the estimable abilities, used to report parameters on a fixed scale. + def scale(abilities, active) + values = abilities.select.with_index { |_, i| active[i] } + return [0.0, 1.0] if values.empty? + + mean, sd = mean_and_sd(values) + [mean, sd.zero? ? 1.0 : sd] + end + def mean_and_sd(values) mean = values.sum / values.size [mean, Math.sqrt(values.sum { |v| (v - mean)**2 } / values.size)] diff --git a/lib/irt_ruby/rasch_model.rb b/lib/irt_ruby/rasch_model.rb index 851c2ec..8f787b5 100644 --- a/lib/irt_ruby/rasch_model.rb +++ b/lib/irt_ruby/rasch_model.rb @@ -3,6 +3,7 @@ require "irt_ruby/estimation" require "irt_ruby/model_options_validator" require "irt_ruby/response_data_validator" +require "irt_ruby/response_graph" module IrtRuby # A class representing the Rasch model for Item Response Theory (ability - difficulty). @@ -11,7 +12,7 @@ module IrtRuby # - Missing data handling (skip nil) # - Multiple convergence checks (log-likelihood + parameter updates) # - Exclusion of extreme response patterns, which have no finite estimates - # - Identification by centering abilities at zero + # - Identification by reporting abilities centered at zero class RaschModel MISSING_STRATEGIES = %i[ignore treat_as_incorrect treat_as_correct].freeze @@ -44,7 +45,7 @@ def initialize(data, @random = seed.nil? ? nil : Random.new(seed) @responses = data_array.map { |row| row.map { |resp| resolve_missing(resp).first } } - @active_persons, @active_items = Estimation.estimable(@responses) + @active_persons, @active_items = ResponseGraph.estimable(@responses) @person_counts, @item_counts = Estimation.response_counts(@responses, @active_persons, @active_items) # Initialize parameters near zero @@ -125,8 +126,8 @@ def average_param_update(old_abilities, old_difficulties) deltas.sum / deltas.size end - # Returns estimates with nil for persons and items excluded as extreme or - # unobserved, plus whether the convergence criteria were met. + # Returns estimates on the identified scale, with nil for persons and items that + # have no finite estimate on it, plus whether the convergence criteria were met. def fit converged = false iterations = 0 @@ -159,12 +160,12 @@ def fit end end - center_abilities end + mean, = Estimation.scale(@abilities, @active_persons) { - abilities: Estimation.estimates(@abilities, @active_persons), - difficulties: Estimation.estimates(@difficulties, @active_items), + abilities: Estimation.estimates(@abilities.map { |theta| theta - mean }, @active_persons), + difficulties: Estimation.estimates(@difficulties.map { |b| b - mean }, @active_items), converged: converged, iterations: iterations } @@ -175,12 +176,5 @@ def stationary? Estimation.stationary?(grad_abilities, [grad_difficulties], @person_counts, @item_counts, @gradient_tolerance) end private :stationary? - - def center_abilities - mean, = Estimation.mean_and_sd(@abilities.select.with_index { |_, i| @active_persons[i] }) - @abilities.map! { |theta| theta - mean } - @difficulties.map! { |b| b - mean } - end - private :center_abilities end end diff --git a/lib/irt_ruby/response_graph.rb b/lib/irt_ruby/response_graph.rb new file mode 100644 index 0000000..d6e8ffb --- /dev/null +++ b/lib/irt_ruby/response_graph.rb @@ -0,0 +1,85 @@ +# frozen_string_literal: true + +module IrtRuby + # Finds the persons and items whose parameters are jointly estimable. + # @api private + module ResponseGraph + module_function + + # Returns [persons, items] flags marking the parameters that have finite joint + # maximum-likelihood estimates on a common scale. In the response graph a correct + # answer links item -> person and an incorrect one person -> item; estimates exist + # within a strongly connected component, so the largest one is kept. All-correct, + # all-incorrect and unobserved rows or columns, perfectly separated subsets, and + # subsets sharing no responses all fall outside it. + def estimable(responses) + num_persons = responses.size + num_items = responses.first.size + forward = Array.new(num_persons + num_items) { [] } + backward = Array.new(num_persons + num_items) { [] } + responses.each_with_index do |row, i| + row.each_with_index do |value, j| + next if value.nil? + + from, to = value == 1 ? [num_persons + j, i] : [i, num_persons + j] + forward[from] << to + backward[to] << from + end + end + + component = strongly_connected_components(forward, backward) + sizes = Hash.new(0) + component.each { |c| sizes[c] += 1 } + largest, size = sizes.max_by { |c, n| [n, -c] } + return [Array.new(num_persons, false), Array.new(num_items, false)] if size < 2 + + [Array.new(num_persons) { |i| component[i] == largest }, + Array.new(num_items) { |j| component[num_persons + j] == largest }] + end + + # Iterative Kosaraju: returns the component index of every node. + def strongly_connected_components(forward, backward) + order = [] + visited = Array.new(forward.size, false) + forward.each_index do |start| + next if visited[start] + + visited[start] = true + stack = [[start, 0]] + until stack.empty? + node, edge = stack.last + if edge < forward[node].size + stack.last[1] += 1 + target = forward[node][edge] + next if visited[target] + + visited[target] = true + stack << [target, 0] + else + order << node + stack.pop + end + end + end + + component = Array.new(forward.size) + count = 0 + order.reverse_each do |start| + next if component[start] + + component[start] = count + stack = [start] + until stack.empty? + backward[stack.pop].each do |source| + next if component[source] + + component[source] = count + stack << source + end + end + count += 1 + end + component + end + end +end diff --git a/lib/irt_ruby/three_parameter_model.rb b/lib/irt_ruby/three_parameter_model.rb index 7fccde1..8c61a77 100644 --- a/lib/irt_ruby/three_parameter_model.rb +++ b/lib/irt_ruby/three_parameter_model.rb @@ -3,6 +3,7 @@ require "irt_ruby/estimation" require "irt_ruby/model_options_validator" require "irt_ruby/response_data_validator" +require "irt_ruby/response_graph" module IrtRuby # A class representing the Three-Parameter model (3PL) for Item Response Theory. @@ -13,7 +14,7 @@ module IrtRuby # - Multiple convergence checks # - Separate gradient calculation & updates # - Exclusion of extreme response patterns, which have no finite estimates - # - Identification by standardizing abilities to mean 0 and SD 1 + # - Identification by reporting abilities standardized to mean 0 and SD 1 class ThreeParameterModel MISSING_STRATEGIES = %i[ignore treat_as_incorrect treat_as_correct].freeze DISCRIMINATION_BOUNDS = [0.01, 5.0].freeze @@ -45,8 +46,9 @@ def initialize(data, @random = seed.nil? ? nil : Random.new(seed) @responses = data_array.map { |row| row.map { |resp| resolve_missing(resp).first } } - @active_persons, @active_items = Estimation.estimable(@responses) + @active_persons, @active_items = ResponseGraph.estimable(@responses) @person_counts, @item_counts = Estimation.response_counts(@responses, @active_persons, @active_items) + @floor_items = floor_items # Initialize parameters @abilities = Array.new(num_rows) { random_between(-0.25..0.25) } @@ -101,6 +103,7 @@ def log_likelihood Math.log(1.0 - c) + Estimation.log_sigmoid(-z) end end + @floor_items.each { |j| ll += Math.log(1.0 - @guessings[j]) } ll end @@ -116,22 +119,40 @@ def compute_gradient b = @difficulties[j] c = @guessings[j] - # With P = c + (1 - c) * s and s = sigmoid(a * (theta - b)): - # dLL/dz = (u - P) * s / P and dLL/dc = (u - P) / ((1 - c) * P). - s = sigmoid(a * (theta - b)) + # With P = c + (1 - c) * s and s = sigmoid(z), z = a * (theta - b): + # dLL/dz = (u - P) * s / P, and dLL/dc = (1 - s) / P if u = 1, -1 / (1 - c) if u = 0. + z = a * (theta - b) + s = sigmoid(z) prob = c + ((1.0 - c) * s) - error = value - prob - dz = prob.zero? ? error : error * s / prob + dz = c.zero? ? value - s : (value - prob) * s / prob grad_abilities[i] += dz * a grad_difficulties[j] -= dz * a grad_discriminations[j] += dz * (theta - b) - grad_guessings[j] += prob.zero? ? 0.0 : error / ((1.0 - c) * prob) + grad_guessings[j] += if value.zero? + -1.0 / (1.0 - c) + else + c.zero? ? Math.exp(-z) : (1.0 - s) / prob + end end + @floor_items.each { |j| grad_guessings[j] -= 1.0 / (1.0 - @guessings[j]) } [grad_abilities, grad_difficulties, grad_discriminations, grad_guessings] end + # Excluded respondents who answered every estimable item incorrectly have ability + # tending to -infinity, where P = c, so each of those responses still contributes + # log(1 - c) to the likelihood of its item's guessing parameter. + def floor_items + @responses.each_index.flat_map do |i| + next [] if @active_persons[i] + + items = @responses[i].each_index.select { |j| @active_items[j] && !@responses[i][j].nil? } + items.all? { |j| @responses[i][j].zero? } ? items : [] + end + end + private :floor_items + def apply_gradient_update(ga, gd, gdisc, gc) old_a = @abilities.dup old_d = @difficulties.dup @@ -165,8 +186,8 @@ def average_param_update(old_a, old_d, old_disc, old_c) deltas.sum / deltas.size end - # Returns estimates with nil for persons and items excluded as extreme or - # unobserved, plus whether the convergence criteria were met. + # Returns estimates on the identified scale, with nil for persons and items that + # have no finite estimate on it, plus whether the convergence criteria were met. def fit converged = false iterations = 0 @@ -200,36 +221,30 @@ def fit end end - standardize_abilities end + mean, sd = Estimation.scale(@abilities, @active_persons) { - abilities: Estimation.estimates(@abilities, @active_persons), - difficulties: Estimation.estimates(@difficulties, @active_items), - discriminations: Estimation.estimates(@discriminations, @active_items), + abilities: Estimation.estimates(@abilities.map { |theta| (theta - mean) / sd }, @active_persons), + difficulties: Estimation.estimates(@difficulties.map { |b| (b - mean) / sd }, @active_items), + discriminations: Estimation.estimates(@discriminations.map { |a| a * sd }, @active_items), guessings: Estimation.estimates(@guessings, @active_items), converged: converged, iterations: iterations } end + # Gradients are compared on the reported scale (theta' = (theta - m) / sd, + # b' = (b - m) / sd, a' = a * sd), where they scale by sd, sd and 1 / sd. def stationary? ga, gd, gdisc, gc = compute_gradient gdisc = Estimation.project(gdisc, @discriminations, *DISCRIMINATION_BOUNDS) gc = Estimation.project(gc, @guessings, *GUESSING_BOUNDS) - Estimation.stationary?(ga, [gd, gdisc, gc], @person_counts, @item_counts, @gradient_tolerance) + _, sd = Estimation.scale(@abilities, @active_persons) + Estimation.stationary?(ga.map { |g| g * sd }, + [gd.map { |g| g * sd }, gdisc.map { |g| g / sd }, gc], + @person_counts, @item_counts, @gradient_tolerance) end private :stationary? - - # Rescales the latent trait without changing the likelihood: theta' = (theta - m) / s, - # b' = (b - m) / s, a' = a * s. - def standardize_abilities - mean, sd = Estimation.mean_and_sd(@abilities.select.with_index { |_, i| @active_persons[i] }) - sd = 1.0 if sd.zero? - @abilities.map! { |theta| (theta - mean) / sd } - @difficulties.map! { |b| (b - mean) / sd } - @discriminations.map! { |a| a * sd } - end - private :standardize_abilities end end diff --git a/lib/irt_ruby/two_parameter_model.rb b/lib/irt_ruby/two_parameter_model.rb index a4058a4..c9ee3bc 100644 --- a/lib/irt_ruby/two_parameter_model.rb +++ b/lib/irt_ruby/two_parameter_model.rb @@ -3,6 +3,7 @@ require "irt_ruby/estimation" require "irt_ruby/model_options_validator" require "irt_ruby/response_data_validator" +require "irt_ruby/response_graph" module IrtRuby # A class representing the Two-Parameter model (2PL) for IRT. @@ -13,7 +14,7 @@ module IrtRuby # - Multiple convergence checks # - Separate gradient calculation & parameter update # - Exclusion of extreme response patterns, which have no finite estimates - # - Identification by standardizing abilities to mean 0 and SD 1 + # - Identification by reporting abilities standardized to mean 0 and SD 1 class TwoParameterModel MISSING_STRATEGIES = %i[ignore treat_as_incorrect treat_as_correct].freeze DISCRIMINATION_BOUNDS = [0.01, 5.0].freeze @@ -39,7 +40,7 @@ def initialize(data, max_iter: 1000, tolerance: 1e-6, param_tolerance: 1e-6, @random = seed.nil? ? nil : Random.new(seed) @responses = data_array.map { |row| row.map { |resp| resolve_missing(resp).first } } - @active_persons, @active_items = Estimation.estimable(@responses) + @active_persons, @active_items = ResponseGraph.estimable(@responses) @person_counts, @item_counts = Estimation.response_counts(@responses, @active_persons, @active_items) # Initialize parameters @@ -131,8 +132,8 @@ def average_param_update(old_a, old_d, old_disc) deltas.sum / deltas.size end - # Returns estimates with nil for persons and items excluded as extreme or - # unobserved, plus whether the convergence criteria were met. + # Returns estimates on the identified scale, with nil for persons and items that + # have no finite estimate on it, plus whether the convergence criteria were met. def fit converged = false iterations = 0 @@ -165,34 +166,28 @@ def fit end end - standardize_abilities end + mean, sd = Estimation.scale(@abilities, @active_persons) { - abilities: Estimation.estimates(@abilities, @active_persons), - difficulties: Estimation.estimates(@difficulties, @active_items), - discriminations: Estimation.estimates(@discriminations, @active_items), + abilities: Estimation.estimates(@abilities.map { |theta| (theta - mean) / sd }, @active_persons), + difficulties: Estimation.estimates(@difficulties.map { |b| (b - mean) / sd }, @active_items), + discriminations: Estimation.estimates(@discriminations.map { |a| a * sd }, @active_items), converged: converged, iterations: iterations } end + # Gradients are compared on the reported scale (theta' = (theta - m) / sd, + # b' = (b - m) / sd, a' = a * sd), where they scale by sd, sd and 1 / sd. def stationary? ga, gd, gdisc = compute_gradient gdisc = Estimation.project(gdisc, @discriminations, *DISCRIMINATION_BOUNDS) - Estimation.stationary?(ga, [gd, gdisc], @person_counts, @item_counts, @gradient_tolerance) + _, sd = Estimation.scale(@abilities, @active_persons) + Estimation.stationary?(ga.map { |g| g * sd }, + [gd.map { |g| g * sd }, gdisc.map { |g| g / sd }], + @person_counts, @item_counts, @gradient_tolerance) end private :stationary? - - # Rescales the latent trait without changing the likelihood: theta' = (theta - m) / s, - # b' = (b - m) / s, a' = a * s. - def standardize_abilities - mean, sd = Estimation.mean_and_sd(@abilities.select.with_index { |_, i| @active_persons[i] }) - sd = 1.0 if sd.zero? - @abilities.map! { |theta| (theta - mean) / sd } - @difficulties.map! { |b| (b - mean) / sd } - @discriminations.map! { |a| a * sd } - end - private :standardize_abilities end end diff --git a/spec/irt_ruby/estimation_spec.rb b/spec/irt_ruby/estimation_spec.rb index 3d1e915..3856098 100644 --- a/spec/irt_ruby/estimation_spec.rb +++ b/spec/irt_ruby/estimation_spec.rb @@ -96,13 +96,53 @@ def simulate(rng, abilities, difficulties, discriminations, guessings) expect(sd).to be_within(1e-12).of(1.0) unless model_class == IrtRuby::RaschModel end - it "keeps the likelihood unchanged when fixing the scale" do - model = model_class.new(mixed_data, max_iter: 1, seed: 1) - model.instance_variable_set(:@max_iter, 0) - before = model.log_likelihood + it "reports rescaled parameters with the same likelihood as the fitted ones" do + model = model_class.new(mixed_data, max_iter: 50, seed: 1) + results = model.fit + fitted = model.log_likelihood + parameters.each { |parameter| model.instance_variable_set("@#{parameter}", results[parameter]) } + + expect(model.log_likelihood).to be_within(1e-9).of(fitted) + end + + it "keeps improving on repeated fit calls" do + model = model_class.new([[1, 0], [0, 1]], seed: 1) + model.fit + first = model.log_likelihood model.fit - expect(model.log_likelihood).to be_within(1e-9).of(before) + expect(model.log_likelihood).to be >= first + expect(model.instance_variable_get(:@learning_rate)).to be > 1e-6 + end + + it "keeps only one of two perfectly separated blocks" do + data = [ + [1, 0, 1, 1], + [0, 1, 1, 1], + [0, 0, 1, 0], + [0, 0, 0, 1] + ] + short = model_class.new(data, max_iter: 200, learning_rate: 0.1, seed: 1).fit + long = model_class.new(data, max_iter: 2000, learning_rate: 0.1, seed: 1).fit + + expect(short[:abilities].count(&:nil?)).to eq(2) + expect(short[:difficulties].count(&:nil?)).to eq(2) + long[:abilities].compact.zip(short[:abilities].compact).each do |long_theta, short_theta| + expect(long_theta).to be_within(1e-3).of(short_theta) + end + end + + it "keeps only one of two blocks that share no responses" do + data = [ + [1, 0, nil, nil], + [0, 1, nil, nil], + [nil, nil, 1, 0], + [nil, nil, 0, 1] + ] + results = model_class.new(data, seed: 1).fit + + expect(results[:abilities].count(&:nil?)).to eq(2) + expect(results[:difficulties].count(&:nil?)).to eq(2) end it "does not report convergence when the learning rate is too small to move" do @@ -141,6 +181,44 @@ def simulate(rng, abilities, difficulties, discriminations, guessings) end end + [IrtRuby::TwoParameterModel, IrtRuby::ThreeParameterModel].each do |model_class| + it "#{model_class} checks stationarity on the reported scale" do + model = model_class.new([[1, 0], [0, 1]], seed: 1) + model.instance_variable_set(:@abilities, [-1e-7, 1e-7]) + model.instance_variable_set(:@difficulties, [0.0, 0.0]) + model.instance_variable_set(:@discriminations, [1.0, 1.0]) + model.instance_variable_set(:@guessings, [0.0, 0.0]) if model_class == IrtRuby::ThreeParameterModel + + expect(model.fit[:converged]).to be(false) + end + end + + describe IrtRuby::ThreeParameterModel do + it "keeps the guessing information of respondents who answered every item incorrectly" do + data = [[1, 0], [0, 1], [0, 0]] + model = described_class.new(data, seed: 1) + model.instance_variable_set(:@abilities, [0.0, 0.0, 0.0]) + model.instance_variable_set(:@difficulties, [0.0, 0.0]) + model.instance_variable_set(:@discriminations, [1.0, 1.0]) + model.instance_variable_set(:@guessings, [0.2, 0.2]) + gradients = model.compute_gradient + + expect(model.log_likelihood).to be_within(1e-12).of(2 * (Math.log(0.6) + Math.log(0.4) + Math.log(0.8))) + expect(gradients[3][0]).to be_within(1e-12).of((0.5 / 0.6) - (2 / 0.8)) + expect(gradients[3][0]).to be_within(1e-6).of(numeric_derivative(model, :guessings, 0)) + end + + it "computes guessing gradients when the item response function underflows" do + model = described_class.new([[1, 0], [0, 1]], seed: 1) + model.instance_variable_set(:@abilities, [-750.0, 0.0]) + model.instance_variable_set(:@difficulties, [0.0, 0.0]) + model.instance_variable_set(:@discriminations, [1.0, 1.0]) + model.instance_variable_set(:@guessings, [0.0, 0.0]) + + expect(model.compute_gradient[3][1]).to be_within(1e-12).of(0.0) + end + end + it "reaches the exact Rasch joint estimates on the README data" do results = IrtRuby::RaschModel.new([[1, 0, 1], [0, 1, 0], [1, 1, 1]], learning_rate: 0.5, seed: 1).fit From 70d9ca40bfffbec4418ea575a23977d4fba69d4c Mon Sep 17 00:00:00 2001 From: Alex Kholodniak Date: Sun, 11 Oct 2026 07:16:13 +0300 Subject: [PATCH 3/6] fix: keep 3PL guessing updates finite and judged on all their responses At c = 0 a badly mispredicted correct response contributes exp(-z) to the guessing gradient, which overflowed to Infinity; backtracking could then never shrink the step, and once the learning rate reached zero the update became NaN and fit raised. Cap each term so the sum stays finite. Guessing gradients also include the limiting terms of excluded all-incorrect respondents, so average them over those responses too when checking stationarity. Document that strong connectivity guarantees finite estimates only for the Rasch model. --- CHANGELOG.md | 2 +- README.md | 2 +- lib/irt_ruby/response_graph.rb | 3 ++- lib/irt_ruby/three_parameter_model.rb | 13 ++++++++++--- spec/irt_ruby/estimation_spec.rb | 21 +++++++++++++++++++++ 5 files changed, 35 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 99985d9..325afce 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,7 +7,7 @@ All notable changes to this project are documented in this file. ### Fixed - Fixed 3PL gradients for ability, difficulty, discrimination, and guessing; they now match finite differences of the log-likelihood. - Only the largest strongly connected part of the response graph is estimated; other respondents and items are returned as `nil` instead of drifting with `max_iter` or keeping random initial values. This covers all-correct, all-incorrect, and unobserved rows or columns as well as perfectly separated or disconnected groups. -- 3PL keeps the limiting `log(1 - c)` contribution of excluded respondents who answered every estimated item incorrectly, and computes guessing gradients stably when the response function underflows. +- 3PL keeps the limiting `log(1 - c)` contribution of excluded respondents who answered every estimated item incorrectly, and keeps guessing gradients finite when the response function underflows. - Log-likelihoods are computed with a numerically stable log-sigmoid instead of adding `1e-15` inside `Math.log`, which flattened the objective for badly mispredicted responses. - Models copy response data on construction and return fresh arrays from `fit`, so later mutation of either does not change model state. diff --git a/README.md b/README.md index 1fe2f14..1c18fbd 100644 --- a/README.md +++ b/README.md @@ -161,7 +161,7 @@ These bounds are model restrictions: a 3PL item with a true lower asymptote abov All models estimate person and item parameters jointly by maximizing the likelihood (joint maximum likelihood, JMLE). -- **Which parameters are estimated.** Think of the responses as a directed graph: a correct answer links the item to the respondent, an incorrect answer links the respondent to the item. Finite estimates on a common scale exist only within a strongly connected part of this graph, so IrtRuby estimates the largest such part and returns `nil` for every other respondent and item. This excludes respondents who answered every observed item correctly or every item incorrectly, items that everyone answered correctly or incorrectly, rows and columns with no observed responses, groups that are perfectly separated from the rest (for example, a group of respondents who answered every item of another group correctly), and groups that share no responses with the rest. Missing-data strategies are applied first, so `:treat_as_correct` and `:treat_as_incorrect` can create such patterns. +- **Which parameters are estimated.** Think of the responses as a directed graph: a correct answer links the item to the respondent, an incorrect answer links the respondent to the item. Finite estimates on a common scale exist only within a strongly connected part of this graph, so IrtRuby estimates the largest such part and returns `nil` for every other respondent and item. This excludes respondents who answered every observed item correctly or every item incorrectly, items that everyone answered correctly or incorrectly, rows and columns with no observed responses, groups that are perfectly separated from the rest (for example, a group of respondents who answered every item of another group correctly), and groups that share no responses with the rest. Missing-data strategies are applied first, so `:treat_as_correct` and `:treat_as_incorrect` can create such patterns. For the Rasch model this condition guarantees finite estimates. For 2PL and 3PL it is necessary but not sufficient: 3PL fits in particular can drift toward boundary solutions in which an item's difficulty keeps growing while its guessing parameter carries the item. Treat very large difficulties, or results with `converged: false`, as unreliable. - **3PL and all-incorrect respondents.** Under 3PL, an excluded respondent who answered every estimated item incorrectly still informs the guessing parameters: as their ability tends to minus infinity each response has probability `1 - c`, and that limiting term stays in the likelihood. - **Scale identification.** Results are reported on a fixed scale: Rasch abilities are centered to mean 0; 2PL and 3PL abilities are standardized to mean 0 and standard deviation 1, with difficulties and discriminations transformed so that response probabilities are unchanged. Only estimated respondents enter these statistics. The model keeps its internal parameters, so calling `fit` again continues from where it stopped; convergence is checked on the reported scale. - **Known limitations.** JMLE item estimates are biased when tests are short, and the bias does not vanish as the number of respondents grows. 2PL and 3PL JMLE is also slow to converge with plain gradient ascent. For high-stakes calibration, compare against marginal maximum likelihood software such as R's `mirt`. diff --git a/lib/irt_ruby/response_graph.rb b/lib/irt_ruby/response_graph.rb index d6e8ffb..b6d50ae 100644 --- a/lib/irt_ruby/response_graph.rb +++ b/lib/irt_ruby/response_graph.rb @@ -11,7 +11,8 @@ module ResponseGraph # answer links item -> person and an incorrect one person -> item; estimates exist # within a strongly connected component, so the largest one is kept. All-correct, # all-incorrect and unobserved rows or columns, perfectly separated subsets, and - # subsets sharing no responses all fall outside it. + # subsets sharing no responses all fall outside it. This suffices for the Rasch + # model; for 2PL and 3PL it is necessary but not sufficient. def estimable(responses) num_persons = responses.size num_items = responses.first.size diff --git a/lib/irt_ruby/three_parameter_model.rb b/lib/irt_ruby/three_parameter_model.rb index 8c61a77..10be79b 100644 --- a/lib/irt_ruby/three_parameter_model.rb +++ b/lib/irt_ruby/three_parameter_model.rb @@ -19,6 +19,10 @@ class ThreeParameterModel MISSING_STRATEGIES = %i[ignore treat_as_incorrect treat_as_correct].freeze DISCRIMINATION_BOUNDS = [0.01, 5.0].freeze GUESSING_BOUNDS = [0.0, 0.35].freeze + # A correct response at c = 0 contributes exp(-z) to the guessing gradient, which + # overflows for very negative z; capping each term keeps the sum finite so that + # backtracking can still shrink the step. + MAX_GUESSING_GRADIENT_TERM = 1e12 def initialize(data, max_iter: 1000, @@ -49,6 +53,8 @@ def initialize(data, @active_persons, @active_items = ResponseGraph.estimable(@responses) @person_counts, @item_counts = Estimation.response_counts(@responses, @active_persons, @active_items) @floor_items = floor_items + @guessing_counts = @item_counts.dup + @floor_items.each { |j| @guessing_counts[j] += 1 } # Initialize parameters @abilities = Array.new(num_rows) { random_between(-0.25..0.25) } @@ -132,7 +138,7 @@ def compute_gradient grad_guessings[j] += if value.zero? -1.0 / (1.0 - c) else - c.zero? ? Math.exp(-z) : (1.0 - s) / prob + c.zero? ? [Math.exp(-z), MAX_GUESSING_GRADIENT_TERM].min : (1.0 - s) / prob end end @floor_items.each { |j| grad_guessings[j] -= 1.0 / (1.0 - @guessings[j]) } @@ -242,8 +248,9 @@ def stationary? gc = Estimation.project(gc, @guessings, *GUESSING_BOUNDS) _, sd = Estimation.scale(@abilities, @active_persons) Estimation.stationary?(ga.map { |g| g * sd }, - [gd.map { |g| g * sd }, gdisc.map { |g| g / sd }, gc], - @person_counts, @item_counts, @gradient_tolerance) + [gd.map { |g| g * sd }, gdisc.map { |g| g / sd }], + @person_counts, @item_counts, @gradient_tolerance) && + Estimation.small_residuals?(gc, @guessing_counts, @gradient_tolerance) end private :stationary? end diff --git a/spec/irt_ruby/estimation_spec.rb b/spec/irt_ruby/estimation_spec.rb index 3856098..b0f20c3 100644 --- a/spec/irt_ruby/estimation_spec.rb +++ b/spec/irt_ruby/estimation_spec.rb @@ -217,6 +217,27 @@ def simulate(rng, abilities, difficulties, discriminations, guessings) expect(model.compute_gradient[3][1]).to be_within(1e-12).of(0.0) end + + it "keeps guessing gradients finite when exp(-z) overflows" do + model = described_class.new([[1, 0], [0, 1]], seed: 1) + model.instance_variable_set(:@abilities, [-750.0, 0.0]) + model.instance_variable_set(:@difficulties, [0.0, 0.0]) + model.instance_variable_set(:@discriminations, [1.0, 1.0]) + model.instance_variable_set(:@guessings, [0.0, 0.0]) + + expect(model.compute_gradient[3]).to all(be_finite) + end + + it "averages guessing gradients over responses of excluded all-incorrect respondents too" do + data = [[1, 0], [0, 1]] + Array.new(100) { [0, 0] } + model = described_class.new(data, learning_rate: 1e-20, seed: 1) + model.instance_variable_set(:@abilities, [-1e-7, 1e-7] + Array.new(100, 0.0)) + model.instance_variable_set(:@difficulties, [472.281575830441] * 2) + model.instance_variable_set(:@discriminations, [0.01] * 2) + model.instance_variable_set(:@guessings, [0.001] * 2) + + expect(model.fit[:converged]).to be(true) + end end it "reaches the exact Rasch joint estimates on the README data" do From e3baa1c098d79e795ae4f3a456d8e05d9d97ef45 Mon Sep 17 00:00:00 2001 From: Alex Kholodniak Date: Sun, 11 Oct 2026 07:19:32 +0300 Subject: [PATCH 4/6] fix: cap guessing gradient terms for tiny positive c too A correct response with a tiny but positive guessing parameter still overflowed through (1 - s) / P, reproducing the NaN crash. Apply the cap to both branches. Strong connectivity does guarantee finite estimates for 2PL with bounded discriminations, so limit the documented caveat to 3PL. --- README.md | 2 +- lib/irt_ruby/response_graph.rb | 3 ++- lib/irt_ruby/three_parameter_model.rb | 8 ++++---- spec/irt_ruby/estimation_spec.rb | 2 ++ 4 files changed, 9 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index 1c18fbd..3a7c7e3 100644 --- a/README.md +++ b/README.md @@ -161,7 +161,7 @@ These bounds are model restrictions: a 3PL item with a true lower asymptote abov All models estimate person and item parameters jointly by maximizing the likelihood (joint maximum likelihood, JMLE). -- **Which parameters are estimated.** Think of the responses as a directed graph: a correct answer links the item to the respondent, an incorrect answer links the respondent to the item. Finite estimates on a common scale exist only within a strongly connected part of this graph, so IrtRuby estimates the largest such part and returns `nil` for every other respondent and item. This excludes respondents who answered every observed item correctly or every item incorrectly, items that everyone answered correctly or incorrectly, rows and columns with no observed responses, groups that are perfectly separated from the rest (for example, a group of respondents who answered every item of another group correctly), and groups that share no responses with the rest. Missing-data strategies are applied first, so `:treat_as_correct` and `:treat_as_incorrect` can create such patterns. For the Rasch model this condition guarantees finite estimates. For 2PL and 3PL it is necessary but not sufficient: 3PL fits in particular can drift toward boundary solutions in which an item's difficulty keeps growing while its guessing parameter carries the item. Treat very large difficulties, or results with `converged: false`, as unreliable. +- **Which parameters are estimated.** Think of the responses as a directed graph: a correct answer links the item to the respondent, an incorrect answer links the respondent to the item. Finite estimates on a common scale exist only within a strongly connected part of this graph, so IrtRuby estimates the largest such part and returns `nil` for every other respondent and item. This excludes respondents who answered every observed item correctly or every item incorrectly, items that everyone answered correctly or incorrectly, rows and columns with no observed responses, groups that are perfectly separated from the rest (for example, a group of respondents who answered every item of another group correctly), and groups that share no responses with the rest. Missing-data strategies are applied first, so `:treat_as_correct` and `:treat_as_incorrect` can create such patterns. For the Rasch model, and for 2PL with its bounded discriminations, this condition guarantees that finite estimates exist (the optimizer may still need more iterations to reach them). For 3PL it is necessary but not sufficient: fits can drift toward boundary solutions in which an item's difficulty keeps growing while its guessing parameter carries the item. Treat very large difficulties, or results with `converged: false`, as unreliable. - **3PL and all-incorrect respondents.** Under 3PL, an excluded respondent who answered every estimated item incorrectly still informs the guessing parameters: as their ability tends to minus infinity each response has probability `1 - c`, and that limiting term stays in the likelihood. - **Scale identification.** Results are reported on a fixed scale: Rasch abilities are centered to mean 0; 2PL and 3PL abilities are standardized to mean 0 and standard deviation 1, with difficulties and discriminations transformed so that response probabilities are unchanged. Only estimated respondents enter these statistics. The model keeps its internal parameters, so calling `fit` again continues from where it stopped; convergence is checked on the reported scale. - **Known limitations.** JMLE item estimates are biased when tests are short, and the bias does not vanish as the number of respondents grows. 2PL and 3PL JMLE is also slow to converge with plain gradient ascent. For high-stakes calibration, compare against marginal maximum likelihood software such as R's `mirt`. diff --git a/lib/irt_ruby/response_graph.rb b/lib/irt_ruby/response_graph.rb index b6d50ae..1de56ab 100644 --- a/lib/irt_ruby/response_graph.rb +++ b/lib/irt_ruby/response_graph.rb @@ -12,7 +12,8 @@ module ResponseGraph # within a strongly connected component, so the largest one is kept. All-correct, # all-incorrect and unobserved rows or columns, perfectly separated subsets, and # subsets sharing no responses all fall outside it. This suffices for the Rasch - # model; for 2PL and 3PL it is necessary but not sufficient. + # model and for 2PL with bounded discriminations; for 3PL it is necessary but not + # sufficient. def estimable(responses) num_persons = responses.size num_items = responses.first.size diff --git a/lib/irt_ruby/three_parameter_model.rb b/lib/irt_ruby/three_parameter_model.rb index 10be79b..dcec658 100644 --- a/lib/irt_ruby/three_parameter_model.rb +++ b/lib/irt_ruby/three_parameter_model.rb @@ -19,9 +19,9 @@ class ThreeParameterModel MISSING_STRATEGIES = %i[ignore treat_as_incorrect treat_as_correct].freeze DISCRIMINATION_BOUNDS = [0.01, 5.0].freeze GUESSING_BOUNDS = [0.0, 0.35].freeze - # A correct response at c = 0 contributes exp(-z) to the guessing gradient, which - # overflows for very negative z; capping each term keeps the sum finite so that - # backtracking can still shrink the step. + # A correct response contributes (1 - s) / P to the guessing gradient (exp(-z) at + # c = 0), which overflows when P is tiny; capping each term keeps the sum finite + # so that backtracking can still shrink the step. MAX_GUESSING_GRADIENT_TERM = 1e12 def initialize(data, @@ -138,7 +138,7 @@ def compute_gradient grad_guessings[j] += if value.zero? -1.0 / (1.0 - c) else - c.zero? ? [Math.exp(-z), MAX_GUESSING_GRADIENT_TERM].min : (1.0 - s) / prob + [c.zero? ? Math.exp(-z) : (1.0 - s) / prob, MAX_GUESSING_GRADIENT_TERM].min end end @floor_items.each { |j| grad_guessings[j] -= 1.0 / (1.0 - @guessings[j]) } diff --git a/spec/irt_ruby/estimation_spec.rb b/spec/irt_ruby/estimation_spec.rb index b0f20c3..6c72350 100644 --- a/spec/irt_ruby/estimation_spec.rb +++ b/spec/irt_ruby/estimation_spec.rb @@ -226,6 +226,8 @@ def simulate(rng, abilities, difficulties, discriminations, guessings) model.instance_variable_set(:@guessings, [0.0, 0.0]) expect(model.compute_gradient[3]).to all(be_finite) + model.instance_variable_set(:@guessings, [1e-310, 0.0]) + expect(model.compute_gradient[3]).to all(be_finite) end it "averages guessing gradients over responses of excluded all-incorrect respondents too" do From b6a0e8f0cf4798ab561473a258b575c5e60d6197 Mon Sep 17 00:00:00 2001 From: Alex Kholodniak Date: Sun, 11 Oct 2026 07:25:01 +0300 Subject: [PATCH 5/6] style: drop matrix alignment spaces flagged by RuboCop 1.92 CI resolves the latest RuboCop because Gemfile.lock is not tracked, and 1.92 reports Layout/ExtraSpacing on these response matrices. --- spec/irt_ruby/three_parameter_model_spec.rb | 12 ++++++------ spec/irt_ruby/two_parameter_model_spec.rb | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/spec/irt_ruby/three_parameter_model_spec.rb b/spec/irt_ruby/three_parameter_model_spec.rb index af704a2..d3bcc6e 100644 --- a/spec/irt_ruby/three_parameter_model_spec.rb +++ b/spec/irt_ruby/three_parameter_model_spec.rb @@ -50,9 +50,9 @@ it "handles nil entries gracefully without raising errors" do missing_data = [ [1, nil, 0], - [1, 0, 1], - [0, 1, nil], - [1, 1, 1] + [1, 0, 1], + [0, 1, nil], + [1, 1, 1] ] model = described_class.new(missing_data, max_iter: 200, learning_rate: 0.05) expect { model.fit }.not_to raise_error @@ -159,9 +159,9 @@ it "handles an entire column missing" do data = [ - [1, nil, 0], - [1, nil, 1], - [0, nil, 1] + [1, nil, 0], + [1, nil, 1], + [0, nil, 1] ] model = described_class.new(data, max_iter: 200) expect { model.fit }.not_to raise_error diff --git a/spec/irt_ruby/two_parameter_model_spec.rb b/spec/irt_ruby/two_parameter_model_spec.rb index 5bebece..0b625e9 100644 --- a/spec/irt_ruby/two_parameter_model_spec.rb +++ b/spec/irt_ruby/two_parameter_model_spec.rb @@ -48,8 +48,8 @@ it "does not raise an error with missing data (nil) in 2PL" do missing_data = [ [1, nil, 0], - [1, 0, 1], - [0, 1, nil] + [1, 0, 1], + [0, 1, nil] ] model = described_class.new(missing_data, max_iter: 200, learning_rate: 0.05) expect { model.fit }.not_to raise_error From e1e1147e5db83ceab3cb0c9e42d79c769cce237d Mon Sep 17 00:00:00 2001 From: Alex Kholodniak Date: Sun, 11 Oct 2026 07:27:24 +0300 Subject: [PATCH 6/6] build: pin RuboCop to 1.92.0 Gemfile.lock is not tracked, so CI resolved whatever RuboCop was newest and, with NewCops enabled, a new release could fail the quality job on unchanged code. RuboCop 1.92 needs Ruby >= 2.7, and Bundler resolves all groups even with BUNDLE_WITHOUT, so the gem is only declared on 2.7+; the Ruby 2.6 job runs just the specs. --- Gemfile | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Gemfile b/Gemfile index ba8e053..c9f03b4 100644 --- a/Gemfile +++ b/Gemfile @@ -8,7 +8,9 @@ gem "rake", "~> 13.0" gem "rspec", "~> 3.0" -gem "rubocop", "~> 1.21" +# Exact pin: NewCops is enabled, so any RuboCop release can add offenses. +# This version needs Ruby >= 2.7; the Ruby 2.6 CI job only runs the specs. +gem "rubocop", "1.92.0", require: false if Gem::Version.new(RUBY_VERSION) >= Gem::Version.new("2.7") # Performance benchmarking gem "benchmark-ips", "~> 2.0"