4
$\begingroup$

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)

$\endgroup$
4
$\begingroup$

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:

enter image description here enter image description here enter image description here

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.

$\endgroup$
  • $\begingroup$ btw I'd like to share this with you if you are interested in PyMC: github.com/CamDavidsonPilon/… $\endgroup$ – Cam.Davidson.Pilon Mar 1 '13 at 16:32
  • $\begingroup$ Thank you very much! I actually have already studied your book, and found it really interesting and understandable. $\endgroup$ – EnricoGiampieri Mar 1 '13 at 16:54
  • $\begingroup$ 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. $\endgroup$ – EnricoGiampieri Mar 1 '13 at 17:00
  • $\begingroup$ 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). $\endgroup$ – Cam.Davidson.Pilon Mar 1 '13 at 17:02
  • 1
    $\begingroup$ 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. $\endgroup$ – EnricoGiampieri Mar 2 '13 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.