8
$\begingroup$

A support case with the identification [CASE:3711872] was created.

EDIT: As of Version 11.1.0 the issue is resolved.


Background

Bayesian Data Analysis, 3rd ed., 2013 by Gelman et al. on page 69 f. has a nice and simple example for using Bayesian inference with a multinomial distribution which has the Dirichlet distribution as its conjugate prior distribution (here the 51st presidential elections in the US). So I felt like checking it quickly.

The Setting

We do find that the posterior distribution for the share of voters supporting Mr. Bush ($x_1$), Mr. Dukakis ($x2$) and what else have you ($x3$) after doing a survey of 1447 adults (observing the proportions $(y1,y2,y3)$) using a noninformative prior distribution (e.g. $Dir(1,1,1)$ ) is given by a Dirichlet distribution:

$ {x1, x2, x3| {y1,y2,y3}} \sim Dir( 728, 584, 138 ) $

In Mathematica this can be encoded as:

posteriorDistr = DirichletDistribution[ { 728, 584, 138 }]

Numerical Issues

We can now do posterior inference about the probability, that Mr. Bush will get more votes than Mr. Dukakis using Probability:

N @ Probability[ x1 > x2, {x1, x2} \[Distributed] posteriorDistr ]

0.999966

Which is the correct result (Mr. Bush is set to become the 41st POTUS). But note what we get, if we use NProbability which is much faster than Probability:

NProbability[ x1 > x2, {x1, x2} \[Distributed] posteriorDistr ]

0.562617

There is no warning and the result is far off the mark.

Why is NProbabilityso bad in choosing an appropriate method in this (still rather simple) case?

(Note, that we do get the correct result, if we give the option Method -> "MonteCarlo" for example.)

EDIT

As of Version 11.1.0 the issue is resolved:

NProbability[x1 > x2, {x1, x2} \[Distributed] posteriorDistr]

0.999964

NProbability will now also pass on a warning that comes up in NIntegrate which tells us that the numerical integration is converging too slowly.

$\endgroup$
  • 6
    $\begingroup$ Just use something like Method -> {"NIntegrate", {MinRecursion -> 2, MaxRecursion -> 5}} in NProbability. $\endgroup$ – ciao Sep 9 '16 at 22:41
7
$\begingroup$

It's more the fault of NIntegrate than NProbability.
Compare

NIntegrate[
 Boole[x1 > x2] PDF[DirichletDistribution[{728, 584, 138}], {x1, x2}], {x1, 0, 1}, {x2, 0, 1}] 

with

N@Integrate[
 Boole[x1 > x2] PDF[DirichletDistribution[{728, 584, 138}], {x1, x2}], {x1, 0, 1}, {x2, 0, 1}]

One can use

NProbability[x1 > x2, {x1, x2} \[Distributed] posteriorDistr, Method -> "Trace"] 

to reveal that NIntegrate is used and the integrand isn't numerically friendly.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Sure thing, but calling NProbability and listening to WRI marketing, you would expect, that there is more internal detection and at least some kind of warning? $\endgroup$ – gwr Sep 10 '16 at 22:15
  • 1
    $\begingroup$ @gwr Anyhow, you should send them a bug report. $\endgroup$ – Karsten 7. Sep 11 '16 at 5:51
  • 2
    $\begingroup$ @Karsten7. +1 but NProbability potentially could analyse the problem and select appropriate options to pass to NIntegrate. It has much more context to work with than NIntegrate $\endgroup$ – mikado Sep 11 '16 at 9:42

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.