# 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!

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. Oct 3, 2019 at 11:21