API for defining parametric probability distributions

Question

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, "Statistics`Library`"]
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

Answer 1

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 *)

Answer 2

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.

Answer 3

$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, "Statistics`Library`"];

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

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

Although fi2 is much slower

fi1 === fi2

(*  True  *)