# Mixture of probits: understanding truncated-based likelihoods

I am trying to implement a mixture model of probits to infer the best decision boundary for every latent subpopulation.

When doing Gibbs sampling, we eventually have to compute $P(y^* | w_c)$ where $y^*$ is the auxilary variable in the augmented version and $w_c$ are the weights or the $\beta$ coefficients in a regression. (for a quick reference to the simple version see Wikipedia. The mixture model does not affect the way $y^*$ and $w$ are sampled.)

Debugging my code, I found that the likelihood of $p(y^* | w_c)$ gives odd results since sometimes they prefer the wrong decision boundary (Test 1). However, the posterior samples pf $w | y^*$ seem to be ok (Test 2)

Test 1 (fail): I sample all the $y^*$ given the true $w_{true}$. Then, given the sampled $y^*$ I compare their likelihood given $w_{true}$, $P(y^* | w_{true})$, with their likelihood given some other $w$, $P(y^* | w_{wrong})$ (in the mixture model, this comparison is done at every gibbs iteration to decide whether this subpopulation should jump to another cluster). What happens is that $P(y^* | w_{wrong}) > P(y^* | w_{true})$.

The original is like this, where the colors represent the value of y* for those with $y>0$ (upper points) and $y<0$ (lower points). The likelihood of this $w$ is -172.0354

Now the likelihood of a clearly wrong $w_c$ is... higher!:

This is another example taken from the real Gibbs sampler. According to this, the third line is more likely than the second one (the true one):

Test 2 (ok): Draws from posterior $w | y^*$. To be sure that given those $y^*$ the most likely $w$ are the ones around the true $w$, I generated samples from the posterior $w | y^*, b, v, y$ where b and v are the prior mean and variance of the normal distribution that generates the w. The result is that no matter my priors b and v, the y* make the posterior to very centered around the true $w$:

My likelihood function (python) is:

logsum = 0
for n in range(len(y_user)):
# if use R function, do not normalize limits a and b
if y_user[n]==1:
low = 0.000001 - dot(X_user[n], w_cluster) # 0 is NOT included in the range
upp = np.inf
else:
low = -np.inf
upp = 0 - dot(X_user[n], w_cluster) # 0 IS included in the range

likelihood = truncnorm.pdf(y_star_user[n], low, upp, loc=dot(X_user[n], w_cluster), scale=1)
logsum += np.log(likelihood)
return logsum


How is it possible that some clearly wrong lines are more likely than the true ones?

Is my likelihood function wrong ?

Update: I normalize the w all the time to avoid numerical issues in the truncated normal function. I think this is a mistake since you avoid the classifier to be very confident about points that are far from the line. But I can't see whether this might be the cause of the problem. Maybe they both the good classifier and the bad are so "shy" that the data is, "by chance" (a missfortunate set of $y^*$) making the bad one win?

• Are you certain you have an intercept term? That could affect the likelihood. This is something I often overlook. – user44764 Jul 16 '14 at 15:57
• I have one which is my first weight $w_0$. I mean, it is drawn from the same multivariate normal distribution than the rest. – alberto Jul 24 '14 at 9:10