[SciPy-User] issue estimating the multivariate Polya distribution: why?
Sat Mar 3 17:51:38 CST 2012
On Sat, Mar 3, 2012 at 6:24 PM, <email@example.com> wrote:
> On Sat, Mar 3, 2012 at 4:05 PM, Emanuele Olivetti
> <firstname.lastname@example.org> wrote:
>> Dear All,
>> I am playing with the multivariate Polya distribution - also
>> known as the Dirichlet compound multinomial distribution. A
>> brief description from wikipedia:
>> I made a straightforward implementation of the probability
>> density function in the log-scale here:
>> together with a straightforward montecarlo estimation (by
>> sampling first from a Dirichlet and then computing the log-likelihood
>> of the multinomial) in the log-scale as well. The log-scale was
>> chosen in order to improve numerical stability.
>> If you run the code liked above you should get these two examples:
>> X: [ 0 50 50]
>> alpha: [ 1 10 10]
>> analytic: -5.22892710577
>> montecarlo -5.23470053651
>> X: [100 0 0]
>> alpha: [ 1 10 10]
>> analytic: -51.737395965
>> montecarlo -93.5266543113
>> As you can see in the first case, i.e. X=[0,50,50], there is excellent
>> agreement between the two implementations while in the second case, i.e.
>> x=[100,0,0], there is a dramatic disagreement. Note that the montecarlo
>> estimate is quite stable and if you change the seed of the random
>> number generator you get numbers not too far from -90.
>> So my question is: where does this issue come from? I cannot see
>> mistakes in the implementation (the code is very simple) and
>> I cannot see the source of numerical instability.
> I don't see anything. I was trying out several different parameters,
> and my only guess is that
> the logaddexp is not precise enough in this case. My results (numpy
> 1.5.1) are even worse.
> The probability that you want to calculate is very low
> For larger values it seems to work fine, but it deteriorates fast when
> the loglikelihood drops below -15 or so (with the versions I have
just an observation
with iterations=1e7 I get much better numbers, which are still way
off. But I don't see why this should matter much, since you are
simulating alpha and not low probability events. (unless lot's of tiny
errors add up in different ways)
> Do you need almost zero probability events?
> In a similar problem with poisson mixtures but without monte carlo, I
> was trying out various ways of rescaling, but didn't come up with
> anything useful.
>> Any hint?
>> SciPy-User mailing list
More information about the SciPy-User