I am running HMC to learn a so called item response model that is common in my practice to estimate student abilities and question difficulties. With only these two groups of parameters, the algorithm seems to work ok, but when I add a third group of parameters, discrimination alphas, the estimation of alphas were not very good. When I add a fourth group of parameters, guessing Cs, the algorithm does not work anymore, showing no learning with acceptance rate 0 from the very beginning of the training. Below you can see the setup of my model and the results (it is a modification of the tutorial of rasch model in the edward library examples).

I am very new to statistical computing and the edward package, so I am not sure if I am setting the model wrong, or the hyperparamters wrong, or it is simply the algorithm not working for this occasion?

“”“Rasch model (Rasch, 1960).”“”

fromfutureimport absolute_import

fromfutureimport division

fromfutureimport print_function

import edward as ed

import matplotlib.pyplot as plt

import numpy as np

import tensorflow as tf

import pandas as pd

from edward.models import Bernoulli, Normal, Empirical,Uniform

from scipy.special import expit

### data

n_students = 300

n_questions = 200

bs = 50000

T = 500

#### define a function to simulate responses from a specified group of items and persons

def build_toy_dataset(n_students, n_questions, n_obs):

student_etas = np.random.normal(size=n_students)

question_etas = np.random.normal(size=n_questions)

question_alphas = np.exp(np.random.normal(loc=0., scale=0.25, size=n_questions))

question_guess = np.random.choice([0,0.2,0.25,0.33,0.5],size = n_questions)student_ids = np.random.choice(range(n_students), n_obs)

question_ids = np.random.choice(range(n_questions), n_obs)logits = question_alphas[question_ids]*(student_etas[student_ids] + question_etas[question_ids])

p = question_guess[question_ids] + (1-question_guess[question_ids])*expit(logits)outcomes = np.random.binomial(1, p, n_obs)

data = pd.DataFrame({‘question_id’: question_ids,

‘student_id’: student_ids,

‘outcomes’: outcomes})return data, student_etas, question_etas, question_alphas, question_guess

#### geneartaion of the true theta and true item difficulty

data, true_s_etas, true_q_etas, true_q_alphas, true_q_guess = build_toy_dataset(

n_students, n_questions, n_obs)

obs = data[‘outcomes’].values

student_ids = data[‘student_id’].values.astype(int)

question_ids = data[‘question_id’].values.astype(int)

### MODEL

#### random effects

student_etas = Normal(loc=0.0, scale= 1.0,

sample_shape=n_students)

question_etas = Normal(loc=0.0, scale= 1.0,

sample_shape=n_questions)log_question_alphas=Normal(loc=0.0, scale= 0.25,

sample_shape=n_questions)question_alphas = tf.exp(log_question_alphas)

question_guess = Uniform(low=0.,high=1.,sample_shape=n_questions)

#### the dependent variable

observation_logodds = tf.multiply(tf.gather(question_alphas, question_ids) ,

tf.gather(student_etas, student_ids) + tf.gather(question_etas, question_ids))

probs_obs = tf.gather(question_guess, question_ids)+(tf.ones(data.shape[0])-tf.gather(question_guess,question_ids))*tf.sigmoid(observation_logodds)outcomes = Bernoulli(probs=probs_obs)

### INFERENCE

qstudents = Empirical(params=tf.get_variable(“q_students/params”,

[T, n_students]))

qquestions = Empirical(params=tf.get_variable(“q_questions/params”,

[T, n_questions]))

qlogalphas = Empirical(params=tf.get_variable(“q_log_alphas/params”,

[T,n_questions]))

qguess = Empirical(params=tf.get_variable(“q_guess/params”,

[T,n_questions]))

latent_vars = {

student_etas: qstudents,

question_etas: qquestions,

log_question_alphas: qlogalphas,

question_guess: qguess

}data = {outcomes: obs}

inference = ed.HMC(latent_vars, data)

inference.run(step_size=0.001)

### report

`f, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4) ax1.set_title('Ability Parameters') ax2.set_title('Difficulty Parameters') ax3.set_title('Log Discrimiation Parameters') ax4.set_title('Guessing Parameters') ax1.set_position([0.0,0,1,1]) ax2.set_position([1.1,0,1,1]) ax3.set_position([2.2,0,1,1]) ax4.set_position([3.3,0,1,1]) ax1.scatter(true_s_etas, qstudents.mean().eval()) ax2.scatter(true_q_etas, qquestions.mean().eval()) ax3.scatter(np.log(true_q_alphas), qlogalphas.mean().eval()) ax4.scatter(true_q_guess, qguess.mean().eval()) plt.show()`