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
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)