Speeding up Probability using Monte Carlo

Question

I have a recipe for generating custom distributions, which I want to use Probability[...] on, but am finding that with more than a few variables it very quickly becomes intractable (it runs for hours). I know that in some circumstances you can approximate these probabilities using a Monte Carlo simulation, but for some reason it is not working in my case---even for very small examples.

For example, with:

dist = 
ProbabilityDistribution[
Piecewise[{{3/4, x1 == 1 && x2 == 0}, {1/4, x1 == 2 && x2 == 1}}],
{x1, 1, 2, 1}, {x2, 0, 1, 1}]

then

Probability[x1 == 1, {x1, x2} \[Distributed] dist]

outputs the answer of ¾ just fine, but

NProbability[x1 == 1, {x1, x2} \[Distributed] dist, 
 Method -> "MonteCarlo"]

outputs the error "NProbability: Unable to generate the necessary samples from [dist]."

Strangely, using an example with only one variable it seems to work. For example, using

dist2 = 
 ProbabilityDistribution[
  Piecewise[{{3/4, x1 == 1}, {1/4, x1 == 2}}],
  {x1, 1, 2, 1}]

then

NProbability[x1 == 1, x1 \[Distributed] dist2, Method -> "MonteCarlo"]

gives an answer that is (approximately) ¾, with no errors.

I know there must be something basic about the "MonteCarlo" method that I'm not understanding, but I've had trouble finding much documentation on it. Any help would be much appreciated!

Answer 1

If you try RandomVariate[dist] you get an error message explaining that there is currently no random sampling method for custom multidimensional distributions. It's an unfortunate limitation.

Since your distribution is basically just a probability table for discrete events, you can work around this limitation by using EmpiricalDistribution:

dist = EmpiricalDistribution[{3/4, 1/4} -> {{1, 0}, {2, 1}}]
NProbability[x1 == 1, {x1, x2} \[Distributed] dist, Method -> "MonteCarlo"]

Note also the necessity to specify both coordinates of the distribution, even though you're only interested in the first coordinate. An alternative is:

NProbability[
  x1 == 1, x1 \[Distributed] MarginalDistribution[dist, 1],
  Method -> "MonteCarlo"
]