# API for defining parametric probability distributions

Mathematica has a lot of machinery for working with predefined probability distributions. It is not clear how to make that machinery work with a new distribution.

Suppose I want to define a brand new distribution

MyDistribution[a, b, c]


What is the minimum I need to specify for all the other machinery to kick in?

Do I just need to specify a PDF or CDF? How do I associate it with MyDistribution?

ProbabilityDistribution does some of what I want, but has limitations. Suppose I define

  MyDistribution[a_, b_, c_] := ProbabilityDistribution[SomeFun[a,b,c,x], {x, -Infinity, Infinity}]


then MyDistribution[a, b, c] will automatically be rewritten to ProbabilityDistribution[SomeFun[a,b,c,x], {x, -Infinity, Infinity}], losing information about what the parameters of the distribution are. This makes it hard to have general purpose Fisher information function that works with user defined distributions.

Here is an example of the information loss I am thinking of. Mathematica has a FisherInformation function that uses the explicit distribution parameters to structure the resulting matrix

AppendTo[$ContextPath, "StatisticsLibrary"] FisherInformation[NormalDistribution[m, s]] // MatrixForm  The row and columns of the resulting matrix correspond to the explicit parameters. It is unclear this function should be extended to work on distributions defined using ProbabilityDistribution • What information about the parameters is being lost? If you type ?MyDistribution, you'll get the complete definition. What are the limitations of ProbabilityDistribution you mention? (I'm not suggesting that there aren't any but rather that it is not obvious what limitations affect what you want to do.) – JimB May 18 '16 at 16:04 • For one thing it loses the canonical order of the parameters. ProbabilityDistribution[1/(E^((m - x)^2/(2*s^2))*(Sqrt[2*Pi]*s)), {x, -Infinity, Infinity}] does not tell you the order of m and s, whereas NormalDistribution[m,s] does. Also SomeFun may additionally refer to d, e, f which are just undefined constants and not actual parameters of the family. After the rewriting this distinction is lost. – Daniel Mahler May 18 '16 at 16:23 • It would really help if you could give a specific example. – JimB May 18 '16 at 16:47 • Well, you are using an undocumented function; with these things, consider yourself lucky if it works, but don't expect too much of them. – J. M.'s ennui May 18 '16 at 18:09 • You should check out the "Create your own distribution" presentation by @Sasha Pavlyk. It is getting a bit old but much of it is still probably relevant wolfram.com/events/technology-conference/2011/… – Andy Ross May 18 '16 at 22:47 ## 3 Answers I'm going to argue that parameter information and paradise are not lost (until there's a specific counter-example). (* Define some function *) someFun[x_, μ_, σ_, a_, b_, c_] := (a/σ) Exp[-(x - μ)^b/(c σ^2)] (* Define a probability density function that depends on someFun and some yet to be given constants *) myDistribution[μ_, σ_] := ProbabilityDistribution[someFun[z, μ, σ, a, b, c], {z, -∞, ∞}, Assumptions -> σ > 0]  Now show the definitions of the functions: ?someFun ?myDistribution  Now work with the probability distribution: (* Set the constants so that myDistribution represents a legitimate probability distribution *) constants = {a -> 1/Sqrt[2 π], b -> 2, c -> 2}; (* Determine some properties of myDistribution *) Mean[myDistribution[3, 4] /. constants] (* 3 *) Variance[myDistribution[3, sigma] /. constants] (* sigma^2 *) Expectation[(r - mean)^3, r \[Distributed] myDistribution[mean, sigma] /. constants] (* 0 *)  • For those too lazy to derive normalization constants: set Method -> "Normalize" in ProbabilityDistribution[]. – J. M.'s ennui May 18 '16 at 17:50 • I have added an example to my question – Daniel Mahler May 18 '16 at 18:00 Not really an answer ... more like an extended comment that is too long for the comment box. But I found the question interesting, for a number of reasons: 1. I did not know that Mma had a FisherInformation function hidden away where you found it - how DID you find it? 2. Your question actually highlights one of my pet dislikes - which is the naming of distributions as black boxes - which is exactly what you are trying to emulate. Even for something as well known as the Normal distribution, there are competing parameterisations. For other distributions, there is not even agreement on what functional form the distribution should take, or any number of competing common forms that co-exist. Even when there is such agreement as to functional form, your own example shows beautifully why using black box names is often just inappropriate. You enter: $$\rightarrow$$  FisherInformation[NormalDistribution[m, s]] // MatrixForm  There is no such thing as the Fisher Information of a Normal distribution. Fisher Information is carried out wrt parameters, not with respect to black boxes. Are the parameters of $$N(\mu, \sigma^2)$$: • $$\mu$$ and $$\sigma$$, or • $$\mu$$ and $$\sigma^2$$? Well, it can be either. And we may want to consider both cases, or something entirely different. Each approach is a different problem, with a different solution. Consider your example of a Normal random variable with pdf $$f(x)$$: (source: tri.org.au) The Fisher Information on $$(\mu, \sigma)$$ is: (source: tri.org.au) The Fisher Information on $$(\mu, \sigma^2)$$ is: (source: tri.org.au) To avoid any confusion, I should note that I am using here the FisherInformation function in the mathStatica add-on to Mathematica, and as declaration, that I am one of the authors. • I found FisherInformation because somebody pointed me to another function in StatisticsLibrary and since I could not find any documentation I just did ?? on the whole library. – Daniel Mahler May 19 '16 at 4:24 • I agree that the way Mathematica does distributions is problematic, but I am trying to work with what I got. – Daniel Mahler May 19 '16 at 4:26 $Version

(*  "10.4.1 for Mac OS X x86 (64-bit) (April 11, 2016)"  *)

dist[m_, s_] =
ProbabilityDistribution[
1/(E^((m - x)^2/(2*s^2))*(Sqrt[2*Pi]*s)), {x, -Infinity, Infinity},
Assumptions -> {s > 0, Element[m, Reals]}];

AppendTo[\$ContextPath, "StatisticsLibrary"];

(fi1 = FisherInformation[NormalDistribution[m, s]]) // MatrixForm


(fi2 = FisherInformation[dist[m, s]]) // MatrixForm


Although fi2 is much slower

fi1 === fi2

(*  True  *)

• Using your example I think I finally see the OP's concern. If dist has the parameters defined in the reverse order, then one obtains the exact same matrix (even though the rows and columns should have been reversed) if the FisherInformation function used the parameters in the same order as the function definition. But it appears that FisherInformation uses a lexicographic order: FisherInformation[dist[m,s]] gives a differently ordered matrix than FisherInformation[dist[z,s]]`. At least that's my current theory subject to change by facts. – JimB May 18 '16 at 18:48