How do I add and subtract probability disributions like real numbers?

Question

I'd like your advice: could you recommend a library that allows you to add/subtract/multiply/divide PDFs (Probability Density Functions) like real numbers?

Behind the scenes, it would have to do a Monte Carlo to work the result out, so I'd probably prefer something fast and efficient, that can take advantage of any GPU in the system.

Update:

This is the sort of C# code I am looking for:

  var a = new Normal(0.0, 1.0); // Creates a PDF with mean=0, std. dev=1.0.
  var b = new Normal(0.0, 2.0); // Creates a PDF with mean=0, std. dev=2.0.
  var x = a + b; // Creates a PDF which is the sum of a and b.
                 // i.e. perform a Monte Carlo by taking thousands of samples 
                 // of a and b to construct the resultant PDF.

Update:

What I'm looking for is a method to implement the algebra on "probability shapes" in The Flaw of Averages by Sam Savage. The video Monte Carlo Simulation in Matlab explains the effect I want - a library to perform math on a series of input distributions.

Update:

Searching for the following will produce info on the appropriate libraries:

• "monte carlo library"
• "monte carlo C++"
• "monte carlo Matlab"
• "monte carlo .NET"

Answer 1

The @Risk Developer Kit allows you to start with a set of probability density functions, then perform algebra on the inputs to get some output, i.e. P = A + B.

The keywords on this page can be used to find other competing offerings, e.g. try searching for:

• "monte carlo simulation model C++"
• "monte carlo simulation model .NET"
• "risk analysis toolkit"
• "distributing fitting capabilties".

Its not all that difficult to code this up in a language such as C++ or .NET. The Monte Carlo portion is probably only about 50 lines of code:

• Read "The Flaw Of Averages" by Sam Savage to understand how you can use algebra on "probability shapes".
• Have some method of generating a "probability shape", either by bootstrapping from some sampled data, or from a pre-determined probability density function, or by using the Math.NET probability library.
• Take 10000 samples from the input probability shapes.
• Do the algebra on the samples, i.e. +, -, /, *, etc, to get 1000 outputs. You can also form a probability tree which implies and, or, etc on the inputs.
• Combine these 10000 outputs into a new "probability shape" by putting the results into 100 discrete "buckets".
• Now that we have a new "probability shape", we can then use that as the input into a new probability tree, or perform an integration to get the area, which converts it back into a hard probability number given some threshold.
• The video Monte Carlo Simulation in Matlab explains this entire process much better than I can.