# Speeding up Probability using Monte Carlo

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!

## 1 Answer

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"
]

• Thanks so much! This is working now. Interestingly, the EmpiricalDistribution method you suggest seems to be ridiculously faster than using ProbabilityDistribution on the discrete pdf I created for my custom distribution. E.g.: a series of computations that took 44,000 seconds to do last night using ProbabilityDistribution took less than 1 second using EmpiricalDistribution. – Kevin Oct 3 '19 at 11:21