Add all sorts of mean and dispersion estimation stuff

 * Delete a lot of comments
 * Remove vestigial dependency on RDataGet
 * Use derivatives stuff from PsychometricsBazaarBase

Author: Frankie Robertson

Project.toml

@@ -28,7 +28,6 @@ PrettyPrinting = "54e16d92-306c-5ea0-a30b-337be88ac337"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 PsychometricsBazaarBase = "b0d9cada-d963-45e9-a4c6-4746243987f1"
 QuickHeaps = "30b38841-0f52-47f8-a5f8-18d5d4064379"
-RDataGet = "a115732e-4334-4ecb-8ea3-f683e7f66d4d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
@@ -68,7 +67,6 @@ PrettyPrinting = "0.4.2"
 PrettyTables = "3"
 PsychometricsBazaarBase = "^0.8.6"
 QuickHeaps = "0.2.2"
-RDataGet = "0.1.0"
 Random = "^1.11"
 Reexport = "1"
 Setfield = "^1"

src/Aggregators/ability_estimator.jl

@@ -26,7 +26,7 @@ function pdf(::LikelihoodAbilityEstimator,
     AbilityLikelihood(tracked_responses)
 end
 
-function power_summary(io::IO, ability_estimator::LikelihoodAbilityEstimator)
+function power_summary(io::IO, ::LikelihoodAbilityEstimator)
     println(io, "Ability likelihood distribution")
 end
 

src/ComputerAdaptiveTesting.jl

@@ -30,6 +30,7 @@ include("./NextItemRules/NextItemRules.jl")
 include("./TerminationConditions.jl")
 
 # Combining / running
+include("./DerivedMeasures.jl")
 include("./Rules.jl")
 include("./Sim/Sim.jl")
 include("./DecisionTree/DecisionTree.jl")
@@ -53,4 +54,4 @@ function require_testext()
     return TestExt
 end
 
-end
+end

src/DerivedMeasures.jl

@@ -0,0 +1,163 @@
+module DerivedMeasures
+
+import PsychometricsBazaarBase: power_summary, GridSummary
+using ..Aggregators: TrackedResponses,
+                     Aggregators,
+                     AbilityIntegrator,
+                     AbilityOptimizer,
+                     AbilityEstimator,
+                     ModeAbilityEstimator,
+                     MeanAbilityEstimator,
+                     LikelihoodAbilityEstimator,
+                     DistributionAbilityEstimator,
+                     get_integrator,
+                     expectation
+using FittedItemBanks: domdims
+using ..NextItemRules: AbilityVariance, compute_criteria, compute_criterion, best_item
+using PsychometricsBazaarBase.Integrators: AnyGridIntegrator, get_grid, normdenom
+using PsychometricsBazaarBase.IndentWrappers: indent
+using PsychometricsBazaarBase: IntegralCoeffs
+using PsychometricsBazaarBase: Differentiation
+
+export PointAndSpreadEstimator, MeanAndStdDevEstimator, LaplaceApproxEstimator, SpreadEstimator
+
+abstract type PointAndSpreadEstimator end
+
+# TODO: These all recalculate everything at the moment, but they should reuse the results generated during the CAT
+
+struct MeanAndStdDevEstimator{
+    DistEstT <: DistributionAbilityEstimator,
+    IntegratorT <: AbilityIntegrator
+} <: PointAndSpreadEstimator
+    dist_est::DistEstT
+    integrator::IntegratorT
+end
+
+MeanAndStdDevEstimator(ability_estimator::MeanAbilityEstimator) = MeanAndStdDevEstimator(ability_estimator.dist_est, ability_estimator.integrator)
+MeanAndStdDevEstimator(ability_variance::AbilityVariance) = MeanAndStdDevEstimator(ability_variance.dist_est, ability_variance.integrator)
+
+function (est::MeanAndStdDevEstimator)(tracked_responses::TrackedResponses)
+    denom = normdenom(est.integrator,
+        est.dist_est,
+        tracked_responses)
+    mean = expectation(IntegralCoeffs.id,
+        domdims(tracked_responses.item_bank),
+        est.integrator,
+        est.dist_est,
+        tracked_responses,
+        denom)
+    return (
+        mean,
+        sqrt(expectation(IntegralCoeffs.SqDev(mean),
+            domdims(tracked_responses.item_bank),
+            est.integrator,
+            est.dist_est,
+            tracked_responses,
+            denom))
+    )
+end
+
+function power_summary(io::IO, est::MeanAndStdDevEstimator)
+    println(io, "Mean and standard deviation estimator")
+    indent_io = indent(io, 2)
+    power_summary(indent_io, est.dist_est)
+    power_summary(indent_io, est.integrator)
+end
+
+show(io::IO, ::MIME"text/plain", est::MeanAndStdDevEstimator) = power_summary(io, est)
+
+struct LaplaceApproxEstimator{
+    DistEstT <: DistributionAbilityEstimator,
+    OptimizerT <: AbilityOptimizer
+} <: PointAndSpreadEstimator
+    dist_est::DistEstT
+    optimizer::OptimizerT
+end
+
+LaplaceApproxEstimator(ability_estimator::ModeAbilityEstimator) = LaplaceApproxEstimator(ability_estimator.dist_est, ability_estimator.optim)
+
+function (est::LaplaceApproxEstimator)(tracked_responses::TrackedResponses)
+    # TODO: Numerical stability: Should directly access the log-pdf here
+    mode = est.optimizer(IntegralCoeffs.one, est.dist_est, tracked_responses)
+    return (
+        mode,
+        -Differentiation.double_derivative((ability -> log(pdf(est, tracked_responses))), mode)
+    )
+end
+
+function power_summary(io::IO, est::LaplaceApproxEstimator)
+    println(io, "Laplace approximation based mean and standard deviation estimator")
+    indent_io = indent(io, 2)
+    power_summary(indent_io, est.dist_est)
+    power_summary(indent_io, est.optimizer)
+end
+
+struct SpreadEstimator{InnerT <: PointAndSpreadEstimator}
+    inner::InnerT
+end
+
+function (est::SpreadEstimator)(tracked_responses::TrackedResponses)
+    _, stddev = est.inner(tracked_responses)
+    return stddev
+end
+
+struct DistributionSampler{
+    DistEst <: DistributionAbilityEstimator,
+    IntegratorT <: AbilityIntegrator,
+    ContainerT <: Union{Vector{Float64}, Vector{Vector{Float64}}}
+}
+    dist_est::DistEst
+    integrator::IntegratorT
+    points::ContainerT
+end
+
+_get_estimator_and_integrator(ability_estimator::MeanAbilityEstimator) = (ability_estimator.dist_est, ability_estimator.integrator)
+_get_estimator_and_integrator(ability_variance::AbilityVariance) = (ability_variance.dist_est, ability_variance.integrator)
+
+function DistributionSampler(composite::Union{MeanAbilityEstimator, AbilityVariance}, points=nothing)
+    dist_est, integrator = _get_estimator_and_integrator(composite)
+    return DistributionSampler(dist_est, integrator, points)
+end
+
+function DistributionSampler(dist_est::DistributionAbilityEstimator, integrator::Union{AbilityIntegrator, Nothing}=nothing, points::Nothing=nothing)
+    @info "DistributionSampler" dist_est integrator points
+    if isnothing(integrator)
+        return nothing
+    end
+    inner_integrator = get_integrator(integrator)
+    if !isnothing(points)
+        return DistributionSampler(dist_est, integrator, points)
+    elseif inner_integrator isa AnyGridIntegrator
+        return DistributionSampler(dist_est, integrator, get_grid(inner_integrator))
+    else
+        return nothing
+    end
+end
+
+function eachmatcol(xs::AbstractMatrix)
+    eachcol(xs)
+end
+
+function eachmatcol(xs::AbstractVector)
+    xs
+end
+
+function (est::DistributionSampler)(tracked_responses::TrackedResponses)
+    num = Aggregators.pdf.(
+        est.dist_est,
+        tracked_responses,
+        eachmatcol(est.points)
+    )
+    denom = normdenom(est.integrator, est.dist_est, tracked_responses)
+    return num ./ denom
+end
+
+function power_summary(io::IO, est::DistributionSampler)
+    println(io, "Distribution sampler")
+    indent_io = indent(io, 2)
+    power_summary(indent_io, est.dist_est)
+    power_summary(indent_io, est.integrator)
+    power_summary(indent_io, GridSummary(est.points))
+end
+
+end

src/NextItemRules/NextItemRules.jl

@@ -22,6 +22,7 @@ using PsychometricsBazaarBase.ConfigTools: @requiresome, @returnsome,
                                            find1_instance, find1_type
 using PsychometricsBazaarBase.Integrators: Integrator, intval
 using PsychometricsBazaarBase: Integrators
+using PsychometricsBazaarBase: Differentiation
 using PsychometricsBazaarBase.IndentWrappers: indent
 import PsychometricsBazaarBase.IntegralCoeffs
 using FittedItemBanks: AbstractItemBank, DiscreteDomain, DomainType,
@@ -59,6 +60,9 @@ export AbilityCovarianceStateMultiCriterion, StateMultiCriterion, ItemMultiCrite
 export InformationMatrixCriteria
 export ScalarizedStateCriterion, ScalarizedItemCriterion
 export DRuleItemCriterion, TRuleItemCriterion
+export ObservedInformationPointwiseItemCriterion
+export RawEmpiricalInformationPointwiseItemCriterion
+export EmpiricalInformationPointwiseItemCriterion
 
 # Prelude
 include("./prelude/abstract.jl")

src/NextItemRules/criteria/pointwise/information.jl

@@ -9,7 +9,7 @@ function compute_criterion(
     ability,
     category
 )
-    actual = -double_derivative((ability -> log_resp(ir, category, ability)), ability) .* resp(ir, category, ability)
+    actual = -Differentiation.double_derivative((ability -> log_resp(ir, category, ability)), ability) .* resp(ir, category, ability)
     -actual
 end
 
@@ -18,7 +18,7 @@ function compute_criterion_vec(
     ir::ItemResponse,
     ability
 )
-    actual = -double_derivative((ability -> log_resp_vec(ir, ability)), ability) .* resp_vec(ir, ability)
+    actual = -Differentiation.double_derivative((ability -> log_resp_vec(ir, ability)), ability) .* resp_vec(ir, ability)
     -actual
 end
 
@@ -87,7 +87,7 @@ function compute_criterion(
         ir,
         ability,
         category
-    ) .- double_derivative((ability -> resp(ir, category, ability)), ability)
+    ) .- Differentiation.double_derivative((ability -> resp(ir, category, ability)), ability)
     -actual
 end
 
@@ -100,7 +100,7 @@ function compute_criterion_vec(
         RawEmpiricalInformationPointwiseItemCategoryCriterion(),
         ir,
         ability
-    ) .- double_derivative((ability -> resp_vec(ir, ability)), ability)
+    ) .- Differentiation.double_derivative((ability -> resp_vec(ir, ability)), ability)
     -actual
 end
 

src/NextItemRules/criteria/pointwise/information_support.jl

@@ -70,18 +70,9 @@ log_resp(ir::ItemResponse{<:GuessAndSlipItemBank}, response, θ) = log(resp(ir,
 log_resp(ir::ItemResponse{<:GuessAndSlipItemBank}, θ) = log(resp(ir, θ))
 log_resp_vec(ir::ItemResponse{<:GuessAndSlipItemBank}, θ) = log.(resp_vec(ir, θ))
 
-function vector_hessian(f, x, n)
-    out = ForwardDiff.jacobian(x -> ForwardDiff.jacobian(f, x), x)
-    return reshape(out, n, n, n)
-end
-
-function double_derivative(f, x)
-    ForwardDiff.derivative(x -> ForwardDiff.derivative(f, x), x)
-end
-
 function expected_item_information(ir::ItemResponse, θ::Number)
     exp_resp = resp_vec(ir, θ)
-    d² = double_derivative((θ -> log_resp_vec(ir, θ)), θ)
+    d² = Differentiation.double_derivative((θ -> log_resp_vec(ir, θ)), θ)
     -sum(exp_resp .* d²)
 end
 
@@ -90,7 +81,7 @@ end
 function expected_item_information(ir::ItemResponse, θ::Vector)
     exp_resp = resp_vec(ir, θ)
     n = domdims(ir.item_bank)
-    hess = vector_hessian(θ -> log_resp_vec(ir, θ), θ, n)
+    hess = Differentiation.vector_hessian(θ -> log_resp_vec(ir, θ), θ, n)
     return -sum(eachslice(hess, dims=1) .* exp_resp)
 end
 
@@ -118,4 +109,4 @@ end
 
 function log_resp(ir::ItemResponse{<:ItemBanks.LogItemBank}, resp, θ)
     log(resp(ItemBanks.inner_ir(ir), resp, θ))
-end
+end

src/NextItemRules/criteria/state/ability_variance.jl

@@ -55,13 +55,13 @@ function compute_criterion(
 )::Float64
     # XXX: Not quite sure about this --- is it useful, the MIRT rules cover this case
     mean = expectation(IntegralCoeffs.id,
-        ndims(tracked_responses.item_bank),
+        domdims(tracked_responses.item_bank),
         criterion.integrator,
         criterion.dist_est,
         tracked_responses,
         denom)
     expectation(IntegralCoeffs.SqDev(mean),
-        ndims(tracked_responses.item_bank),
+        domdims(tracked_responses.item_bank),
         criterion.integrator,
         criterion.dist_est,
         tracked_responses,

src/NextItemRules/porcelain/porcelain.jl

@@ -9,3 +9,15 @@ function TRuleItemCriterion(ability_estimator)
         InformationMatrixCriteria(ability_estimator),
         TraceScalarizer())
 end
+
+function ObservedInformationPointwiseItemCriterion()
+    TotalItemInformation(ObservedInformationPointwiseItemCategoryCriterion())
+end
+
+function RawEmpiricalInformationPointwiseItemCriterion()
+    TotalItemInformation(RawEmpiricalInformationPointwiseItemCategoryCriterion())
+end
+
+function EmpiricalInformationPointwiseItemCriterion()
+    TotalItemInformation(EmpiricalInformationPointwiseItemCategoryCriterion())
+end

src/Sim/Sim.jl

@@ -7,6 +7,7 @@ using DocStringExtensions
 using StatsBase
 using FittedItemBanks: AbstractItemBank, ResponseType, ItemResponse, domdims
 using PsychometricsBazaarBase: show_into_buf, power_summary_into_buf
+using PsychometricsBazaarBase.ConfigTools: @requiresome, @returnsome
 using PsychometricsBazaarBase.Integrators
 using PsychometricsBazaarBase.IndentWrappers: indent
 using PsychometricsBazaarBase: GridSummary
@@ -25,6 +26,7 @@ using ..Aggregators: TrackedResponses,
                      MeanAbilityEstimator,
                      LikelihoodAbilityEstimator,
                      RiemannEnumerationIntegrator
+using ..DerivedMeasures: MeanAndStdDevEstimator, LaplaceApproxEstimator, DistributionSampler
 using ..NextItemRules: AbilityVariance, compute_criteria, compute_criterion, best_item
 import Base: show
 import PsychometricsBazaarBase: power_summary