# Bernoulli variable on pymc

Im not fully sure that this is the right place to ask, but I have a problem with pymc that I'm not able to grasp.

I'm trying to simulate a simple counting under two different scenario: Under the first hypothesis I expect a count of ten and under the alternative a count of 30. Given that I presume a probability for the first hypothesis of ~10%, I would like to know the new value after an observation of 25.

The code that I tried to write is this one

frequency = pymc.Beta('frequency', alpha=3.5 , beta=31.5, value=0.5)
gene_variant = pymc.Bernoulli('gene_variant', p=frequency, value=array([False]))

mean = pymc.Lambda('mean', lambda gene_variant=gene_variant: where(gene_variant, 10.0, 30.0))

counting = pymc.Poisson('counting', mu=mean, value=25), observed=True)

model = pymc.MCMC([frequency, gene_variant, mean, counting])
model.sample(iter=4000, burn=500)


A similar approach works fine when trying to model a mixture of distribution, but when I try on this simple problem the value of gene_variant is not updated. When I run the pymc.Matplot.plot diagnostic plot I see a constant value of gene_variant (and thus of mean), no matter how long I let the simulation run.

I am quite confused, as trying the random function of gene_variant I obtain a nice chain of positive and negative results.

It's a very simple problem, so I can't believe it's a bug of pymc, but I can't spot the error. Someone can explain me what's I'm doing wrong?

I am using pymc 2.2 and numpy 1.7

## EDIT:

### model specification

In a population there are two variant of a gene, divided in the population in a 90% / 10% division. This frequency is known with a precision of +- 5%. We know that the result of a biological exam of counting the amount of white cells in the blood for the wild type gene (the most frequent) has a mean of 10, while the variant has a mean of 30. What I want to know is the estimated probability that a person that has a count of 25 is a variant type gene.

To model this I choose the frequency of the gene in the population to be modeled with a Beta distribution with a and b calculated to reflect the know distribution (frequency). The genotype of the person under study is represented by a Bernoulli variable driven by the population frequency (*gene_variant*).

Depending on the genetype we expect a different mean of the count, (mean) that we expect to be modeled by a Poisson variable (counting).

Given only one observation we don't expect to see a variation on the population frequency, but to observe increase of knowledge about the genotype of the person under study (wild type Vs variant)

I ran your code, and performed some diagnostics/debugging. Here's is what I generated, and my interpretation. For transparency, here's my code

frequency = pymc.Beta('frequency', alpha=3.5 , beta=31.5, value=0.5)

gene_variant = pymc.Bernoulli('gene_variant', p=frequency, value=array([False]))

#mean = pymc.Lambda('mean', lambda gene_variant=gene_variant: where(gene_variant, 10.0, 30.0))
"""I prefer to use the boiler plate deterministic template. Plus where confuses me sometimes. Please confirm this is indeed what you intended in your Lambda.
@pymc.deterministic
def mean( gene_variant = gene_variant ):
if gene_variant:
return 10.
else:
return 30.

counting = pymc.Poisson('counting', mu=mean, value=25, observed=True)

model = pymc.MCMC([frequency, gene_variant, mean, counting])
model.sample(iter=10000, burn=5000)

mcplot(model)


Below are my results, let me know if they match yours, more or less:

## Interpretation

What is observed in the above plots is not that unlikely, given the prior and observation. Considering that you assigned a Beta distribtion skewed heavily towards the true mean being 30 and you observe a value of 25 (which under H1 has probability five orders of magnitude less than H2), I would say the results look right: It would be VERY unlikely to have the mean by 10, hence the trace of mean will VERY unlikely be anything but 30 (and consequently gene_variant would be always False.

I did the same analysis with a less strict prior on the frequency ( I chose Uniform over 0,1), and did observe some occurrences of gene_variant == True.

• btw I'd like to share this with you if you are interested in PyMC: github.com/CamDavidsonPilon/… – Cam.Davidson.Pilon Mar 1 '13 at 16:32
• Thank you very much! I actually have already studied your book, and found it really interesting and understandable. – EnricoGiampieri Mar 1 '13 at 16:54
• I guess that my problem (after some trials with the same mean on both distribution) is that the MonteCarlo algorithm is quite slow in the exporation of the parameter space when a bernoulli is involved. May I ask you how to change the algorithm? by the way yes, the lambda and the function should be equivalent. – EnricoGiampieri Mar 1 '13 at 17:00
• Change the algorithm to promote searching the parameter space? I would suggest adding a less restrictive prior, like `frequency = pymc.Uniform('frequency', 0, 1)'. Other wise there is not much one can do: the observations really drive the inference here (see my point on comparing probabilities under H1 vs H2). – Cam.Davidson.Pilon Mar 1 '13 at 17:02
• it would be a pleasure. I will edit the question with a detailed explanation later this afternoon, but in few words is a model of diagnosis of a gene variant by a medical test. – EnricoGiampieri Mar 2 '13 at 15:24