Hi everybody. I read silently the discourse for a few months, and now I’ve decided to write, because I need help with a project.
I am trying to setup a hierarchical model to infer bias and variance for a community of reviewers.
In order to be sure to build the model correctly, I’ve started creating simulated data, asking the model to infer the (known) parameters that generated such data.
- I start with a matrix 𝐸 of size 𝑅×𝐽 with 𝑅 the number of reviewers and 𝐽 the number of reviewed jobs.
- The element 𝐸𝑖𝑗 of matrix contains is a continous number 𝑛 (𝑛≥0) representing the score given by the reviewer 𝑖 to the job 𝑗
- So far every job is judged by just 1 of the 𝑅 reviewers, but for the future I need the model to accept more reviews for a single job
- I’ve created data with 𝑅=5 and 𝐽=1000
Following this implementation of Dawid-Skene model, I created 3 vectors, with 𝐽 elements each:
- train_jobs: contains the idxs of the reviewed jobs
- train_revs: contains the idxs of the reviewers that judged the respective job
- train_evals: contains the score given by the reviewers to the jobs
So that, for each row element k, I have train_revs[k] = 𝑖, train_jobs[k] = 𝑗, train_evals[k] = 𝐸𝑖𝑗.
This is my model:
with pm.Model() as rev_model: alpha_job = pm.Gamma('alpha_job', alpha=1., beta=1./3., shape=J) beta_rev = pm.Uniform('beta_rev', lower=-10, upper=10, shape=R) sigma_rev = pm.HalfNormal('sigma_rev', sigma=1, shape=R) mu_obs = alpha_job[train_jobs] + beta_rev[train_revs] sigma_obs = sigma_rev[train_revs] obs = pm.Normal('epts', mu=mu_obs, sigma=sigma_obs, observed=train_evals) mcmc_trace = pm.sample(chains=4, draws=1000)
I’ve built the dataset so that I have about 200 jobs per reviewer, to make sure that the model could see data for each reviewer.
When I run the sampling I get the errors:
There were 839 divergences after tuning. Increase
The acceptance probability does not match the target. It is 0.2170325472246399, but should be close to 0.8. Try to increase the number of tuning steps.
There were 8 divergences after tuning. Increase
The acceptance probability does not match the target. It is 0.9237895744308975, but should be close to 0.8. Try to increase the number of tuning steps.
There were 15 divergences after tuning. Increase
The acceptance probability does not match the target. It is 0.9145047463020813, but should be close to 0.8. Try to increase the number of tuning steps.
There were 237 divergences after tuning. Increase
The acceptance probability does not match the target. It is 0.5561411316516693, but should be close to 0.8. Try to increase the number of tuning steps.
The rhat statistic is larger than 1.4 for some parameters. The sampler did not converge.
The estimated number of effective samples is smaller than 200 for some parameters.
and I get these parameters
(here I’ve added the first column to the summary df, in which you can read the real value of the “beta_rev” parameters)
I cannot understand what is wrong and why the model seems to be able to infer some parameters and not the others.
To complete my explanation I’ll paste also the code that I use to create the simulated data:
def simulate_data(N_reviewers=10, N_jobs=100, rev_per_job=3, seed=14, sigma_noise=1.): if seed != None: np.random.seed(seed) # sistematic error of reviewers R_biases = np.array([np.random.uniform(-5, 5, N_reviewers)]).T # random error of J_real = np.random.gamma(shape=1, scale=3., size=N_jobs) #J_real = np.random.lognormal(sigma=1.0, size=N_jobs) # biased judgement of jobs EPTs EPTS = R_biases + J_real # keep only rev_per_job reviews for each job for j in range(N_jobs): reviews = np.random.choice(range(N_reviewers), replace=False, size=rev_per_job) remove = set(range(N_reviewers)) - set(reviews) EPTS[list(remove), j] = np.nan # biased and noised values of jobs EPTs EPTS += np.random.normal(loc=0., scale=sigma_noise, size=EPTS.shape) # remove negative EPTs EPTS[np.where(EPTS < 0)] = 0 return R_biases, J_real, EPTS # create E matrix (EPTS) and beta_rev values (biases) N_reviewers=5 N_jobs=1000 rev_per_job=2 R_biases, J_real, EPTS = simulate_data(rev_per_job=rev_per_job, N_reviewers=N_reviewers, N_jobs=N_jobs, seed=14, sigma_noise=1.) # create mask to make the model see just one review np.random.seed(14) mask = np.array([False]*len(all_reviewers)) for i in range(0, len(all_reviewers), 2): rnd = np.random.choice([0, 1]) mask[i + rnd] = True # create the 3 vectors train_revs = all_reviewers[mask] train_jobs = all_jobs[mask] train_evals = all_evaluations[mask] test_revs = all_reviewers[~mask] test_jobs = all_jobs[~mask] test_evals = all_evaluations[~mask]
Here I’ve created the matrix 𝐸 with 2 reviews per job, but then I’ve masked one of the 2 reviews.
Thanks for your help!