0

I would like to estimate IRT model using PyMC3. I generated data with the following distribution:

alpha_fix = 4
beta_fix = 100
theta= np.random.normal(100,15,1000)
prob = np.exp(alpha_fix*(theta-beta_fix))/(1+np.exp(alpha_fix*(theta-beta_fix)))
prob_tt = tt._shared(prob)

Then I created a model using PyMC3 to infer the parameter:

irt = pm.Model()
with irt:
#     Priors
alpha = pm.Normal('alpha',mu = 4 , tau = 1)
beta = pm.Normal('beta',mu = 100 , tau = 15)
thau = pm.Normal('thau' ,mu = 100 , tau = 15)

#     Modelling
p = pm.Deterministic('p',tt.exp(alpha*(thau-beta))/(1+tt.exp(alpha*(thau-beta))))

out = pm.Normal('o',p,observed = prob_tt)

Then I infer through the model:

with irt: 
    mean_field = pm.fit(10000,method='advi', callbacks=[pm.callbacks.CheckParametersConvergence(diff='absolute')])

Finally, Sample from the model to get compute posterior:

pm.plot_posterior(mean_field.sample(1000), color='LightSeaGreen');

But the results of the "alpha" (mean of 2.2) is relatively far from the expected one (4) even though the prior on alpha was well-calibrated.

Posterior parameters

Would you have an idea of the origin of this gap and how to fix it?

Thanks a lot,

Pbcacao
  • 35
  • 4
  • What is the `tau` or `sd` for the measurement Normal (`o`)? Could misspecification of the measurement error could have an effect on the `alpha`? I.e., lack of identifiability. – merv Apr 26 '20 at 00:43
  • I think you need to investigate the convergence of the variational inference as it probably hasn't converged. If you don't have to use VI then this model will recover the value of ```alpha_fix``` by sampling using NUTS: just add the line ```trace = pm.sample()``` after the definition of ```out```. – LeoC Apr 27 '20 at 10:33

1 Answers1

0

out = pm.Normal('o',p,observed = prob_tt)

Why you are using Normal instead of Bernoulli ? Also, what is the variance of normal ?

piapple
  • 1
  • 1