# Looping through random variables
**URL:** https://discourse.pymc.io/t/looping-through-random-variables/2555
**Category:** Questions
**Created:** [January 20, 2019, 8:59am UTC](https://discourse.pymc.io/t/looping-through-random-variables/2555 "2019-01-20T08:59:50Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@Paras\_Chopra](https://discourse.pymc.io/u/Paras_Chopra)
#### Post date: [January 20, 2019, 8:59am UTC](https://discourse.pymc.io/t/looping-through-random-variables/2555/1 "2019-01-20T08:59:50Z")
Hello,
I’m trying to follow Chapter 2 of Model Based Machine Learning [Model-Based Machine Learning: Chapter 2. Assessing People's Skills](http://mbmlbook.com/LearningSkills_Moving_to_real_data.html)
It’s about assessing people’s skills through questions they answer. The generative model is such that there are 7 skills and 48 questions. Each question may test one or more skills. I initialize each skill as a Bernoulli (shape of 7, since there are 7 skills)
> skills = pm.Bernoulli(‘skills’, p=0.5\*np.ones(len(df\_skills.columns)), shape=len(df\_skills.columns))
Next, I want to loop through each question, checking which all skills are required to answer that question. I do that through
```
hasAllSkillsForQuestionP = np.ones(len(df_correct.columns))
for q in range(len(df_correct.columns)):
skills_indexes_needed_for_this_q = np.where(df_skills_per_q.transpose()[q].tolist())[0]
for skill in skills_indexes_needed_for_this_q:
hasAllSkillsForQuestionP[q] = hasAllSkillsForQuestionP[q]*skills[skill]
```
And then I add unpredictability (that is, if a candidate has all skills for answering a question, probability of correct=0.9, if the candidate doesn’t have all skills, probability of correct = guessing = 0.2… it’s MCQ of 5 questions)
> ```
> isCorrectQuestion = pm.Deterministic('isCorrectQuestion', pm.math.switch(hasAllSkillsForQuestionP, pm.Bernoulli(p=0.9), pm.Bernoulli(p=0.2)))
>
> ```
When I try compiling the code, I get this error:
> * * *
ValueError Traceback (most recent call last)
in
15 for skill in skills\_indexes\_needed\_for\_this\_q:
16
—\> 17 hasAllSkillsForQuestionP[q] = hasAllSkillsForQuestionP[q]\*skills[skill]
18
19 #hasAllSkillsForQuestion = pm.Bernoulli(‘hasAllSkillsForQuestion’, p=hasAllSkillsForQuestionP, shape=len(hasAllSkillsForQuestionP))
**ValueError: setting an array element with a sequence.**
Why is this happening? As per [docs](https://docs.pymc.io/notebooks/api_quickstart.html), random variables with shape support full indexing like
> with model:
> y = x[0] \* x[1] # full indexing is supported
How do I fix this?
---
### Author:  [@junpenglao](https://discourse.pymc.io/u/junpenglao)
#### Post date: [January 20, 2019, 5:35pm UTC](https://discourse.pymc.io/t/looping-through-random-variables/2555/2 "2019-01-20T17:35:11Z")
You cannot assign theano tensor to `hasAllSkillsForQuestionP` which is a numpy array. You can either create an empty theano tensor and [assign value to it](https://stackoverflow.com/questions/32280071/how-to-assign-values-elementwise-to-theano-matrix-difference-between-numpy-and), or index to `skills` and reshape it.
---
### Author:  [@Paras\_Chopra](https://discourse.pymc.io/u/Paras_Chopra)
#### Post date: [January 21, 2019, 7:23am UTC](https://discourse.pymc.io/t/looping-through-random-variables/2555/3 "2019-01-21T07:23:53Z")
Thanks. I changed
> hasAllSkillsForQuestionP = np.ones(len(df\_correct.columns))
to
> hasAllSkillsForQuestionP = pm.Deterministic(‘hasAllSkillsForQuestionP’,theano.shared(np.ones((22, len(df\_correct.columns)))))
and
> hasAllSkillsForQuestionP[q] = hasAllSkillsForQuestionP[q]\*skills[skill]
to
> tt.set\_subtensor(hasAllSkillsForQuestionP[:,q], hasAllSkillsForQuestionP[:,q]\*skills[:,int(skill)])
where (tt is theano library) and additional dimension is what I added later for my model (not related to my question).
The model compiled successfully. Thanks for your help!
---
### Author:  [@junpenglao](https://discourse.pymc.io/u/junpenglao)
#### Post date: [January 21, 2019, 9:31am UTC](https://discourse.pymc.io/t/looping-through-random-variables/2555/4 "2019-01-21T09:31:45Z")
Glad that you get it works. Still you might want to avoid using for loop to improve performance.
---
### Author:  [@Paras\_Chopra](https://discourse.pymc.io/u/Paras_Chopra)
#### Post date: [January 21, 2019, 10:31am UTC](https://discourse.pymc.io/t/looping-through-random-variables/2555/5 "2019-01-21T10:31:47Z")
Thanks. Yes, that’s what I ended up doing. Instead of looping, I did matrix multiplication (which had the same effect).
My updated model looks like this:
```
skills = pm.Bernoulli('skills', p=0.5*np.ones((22, 7)), shape=(22,7))
hasAllSkillsForQuestionP = tt.dot(skills, df_skills_per_q_t.T.values)
hasAllSkillsForQuestion = pm.Bernoulli('hasAllSkillsForQuestion', p=hasAllSkillsForQuestionP, shape=(22,len(df_correct.columns)))
isCorrectQuestion = pm.Bernoulli('isCorrectQuestion', p=pm.math.switch(hasAllSkillsForQuestion, 0.9, 0.2), observed=df_correct, shape=(22, 48))
```
But now the problem has shifted to inference. When I run inference for 1000 steps, it selects BinaryGibbsMetropolis: [skills, hasAllSkillsForQuestion] and after 1000 steps, when I inspect the trace, it’s really bad. Most posteriors have NaN as n\_eff (which I interpret as effective samples to predict the value). Although some have 3 or 4 n\_eff.
Do you know what might be happening, and how do I do better inference?
In the book they use message passing algorithms for inference, but I understand that PyMC3 doesn’t implement it.
---
### Author:  [@junpenglao](https://discourse.pymc.io/u/junpenglao)
#### Post date: [January 21, 2019, 10:38am UTC](https://discourse.pymc.io/t/looping-through-random-variables/2555/6 "2019-01-21T10:38:57Z")
In your model you have a latent discrete variable, which is why sampling needs a `BinaryGibbsMetropolis`. In general, discrete latent variable are very difficult to inference, the general recommendation is that to treat it as a mixture model problem and rewrite it into a marginalized mixture model.
I think this post is the closest to your problem that you can start with:
> [@Problem with hierachical occupancy model](https://discourse.pymc.io/t/problem-with-hierachical-occupancy-model/1702/6):
>
> I would do it this way: First reshape everything so within the pm.Model we are not dealing with dict: index = [] W\_r = [] y\_r = [] for i in range(n): index.append(np.ones((W[i].shape[0], 1))\*i) W\_r.append(W[i]) y\_r.append(y[i][:, np.newaxis]) index = np.concatenate(index, axis=0).astype(int64) W\_r = np.concatenate(W\_r, axis=0) y\_r = np.concatenate(y\_r, axis=0) Then I will use matrix multiplication instead of nested for loops: with pm.Model() as m: beta = pm.Normal('beta', mu=…
See also:
> [@Naive Bayes model with PyMC3](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/10):
>
> So here is how I would approach it: First, I would rewrite the explicit model, so that the observed is a flatten array - it is easier to handle shape wise: # flatten and index data data\_flatten = np.reshape(data, data.shape[0]\*data.shape[1]) data\_index = np.repeat(np.arange(data.shape[0]), data.shape[1]) with pm.Model() as model: # Global topic distribution theta = pm.Dirichlet("theta", a=alpha) # Word distributions for K topics phi = pm.Dirichlet("phi", a=beta, shape=(K,…
It provides a way to marginalizing the latent discrete variables.
---
### Author:  [@Paras\_Chopra](https://discourse.pymc.io/u/Paras_Chopra)
#### Post date: [January 21, 2019, 1:45pm UTC](https://discourse.pymc.io/t/looping-through-random-variables/2555/7 "2019-01-21T13:45:37Z")
Thanks. I do not fully understand mixture models or how to construct marginalized models. Is there a tutorial that I should follow?
By the way, I end up changing Bernoulli
> skills = pm.Bernoulli(‘skills’, p=0.5\*np.ones((22, 7)), shape=(22,7))
To Beta distribution
> skills = pm.Beta(‘skills’, alpha=2, beta=2, shape=(22,7))
And model is providing better estimates. I also felt perhaps I should have used beta in the first place as probability of skill could be between 0 or 1, rather than being binary.