# Given an exact formula, how can Mathematica find a probability distribution whose PDF matches it?

So, given some data, Mathematica 10.2 can now attempt to figure out what probability distribution might have produced it. Cool! But suppose that, instead of having data, we have something that is in some ways better -- a formula. Let's call it $f$. We suspect -- perhaps because $f$ is non-negative over some domain and because the integral of $f$ over that domain is 1 -- that $f$ is actually the PDF of some distribution (Normal, Lognormal, Gamma, Weibull, etc.) or some relatively simple transform of that distribution.

Is there any way that Mathematica can help figure out the distribution (or simple transform) whose PDF is the same as $f$?

Example: Consider the following formula:

1/(2*E^((-m + Log[5])^2/8)*Sqrt[2*Pi])


$$\frac{e^{-\frac{1}{8} (\log (5)-m)^2}}{2 \sqrt{2 \pi }}$$

As it happens -- and as I discovered with some research and guesswork -- this formula is the PDF of NormalDistribution[Log[5], 2] evaluated at $m$. But is there a better way than staring or guessing to discover this fact? That is, help me write FindExactDistribution[f_, params_].

## Notes

• The motivation for the problem comes from thinking about Conjugate Prior distributions but I suspect it might have a more general application.

• One could start with mapping PDF evaluated at $m$ over a variety of continuous distributions. And if I did this I would at some point get to what I will call $g$, which is the PDF or the NormalDistribution with parameters $a$ and $b$ evaluated at $m$.

1/(b*E^((-a + m)^2/(2*b^2))*Sqrt[2*Pi])

$$\frac{e^{-\frac{(m-a)^2}{2 b^2}}}{\sqrt{2 \pi } b}$$

But unless I knew that if I replaced $a$ by Log[5] and $b$ by $2$ that I would get $f$, this fact would not mean a lot to me. I suppose I could look at the TreeForm of $f$ and $g$ and I would notice certain similarities, and that might be a hint, but I am not sure how to make much progress beyond that observation. Ultimately, the problem looks to be about finding substitutions in parts of a tree ($g$) which, after evaluation, yield a tree that matches a target $f$. I have the suspicion that this is a difficult problem with an NKS flavor but one for which Mathematica and its ability to transform expressions might be well suited.

I appreciate the responses here. But let me provide an example that is perhaps not so easy. Suppose the target function f is as follows: $\frac{7}{10 (a-2)^2}$ for the domain ($-\infty,\frac{13}{10}$]. If we create a probability distribution out of this and then generate 10,000 random samples from the distribution and then run FindDistribution

 dis = ProbabilityDistribution[7/(10 (-2 + a)^2), {a, -\[Infinity], 13/10}];
rv = RandomVariate[dis,10^4];
fd=FindDistribution[rv,5]


The result is a mixture distribution of normal distributions, a beta distribution, a weibull distribution, a normal distribution and a mixture distribution of a normal distribution and a gamma distribution.

The mixture distributions are clearly of the wrong form, the normal distribution is clearly not right, Although I am not positive, I don't believe the Weibull Distribution or the Beta Distribution is correct either. In fact, I don't know what the correct answer is, though I think it might be a fairly simple transform of a single parameter distribution. The point, however, is that the FindDistribution process, does not seem to work in this case. And that's why I am hoping for something better.

• What if you sampled (with a random number generator) from f assuming that it is a pdf, then used the resulting data in Mathematica to see if it will find a pdf that matches the data? Aug 11, 2015 at 23:26
• Given that the software seems poised to implement almost every distribution there is, one would hope there is a way to find the closest-match among built-ins… Aug 12, 2015 at 3:28
• > "In fact, I don't know what the correct answer is," .... ///////////// I like your new example. It is the pdf of $Y = 2-X$ where $X \sim Pareto(1, 7/10)$. Not sure it would be easy to do that automagically, however. Aug 12, 2015 at 17:07
• There are infinitely many distributions ... (and thus also infinitely many possible matches) Aug 12, 2015 at 17:26
• ... which really highlights the difficulty of providing a general solution to such problems. The number of named special distributions is tiny compared to the number of common functional forms .. and even if one just wishes to restrict focus on named cases, how does one deal with $X^2$ or ratios of common forms, or truncated or censored functions of common forms, or shifted cases, or ... Aug 12, 2015 at 19:05

UPDATE: quite interesting parallel discussion and solutions (see Emerson Willard answer) can be found HERE.

Maybe this is not exactly what you are looking for, but at least this gives you a very close guess and it is easy to figure out the rest.

dis = ProbabilityDistribution[
1/(2*E^((-m + Log[5])^2/8)*Sqrt[2*Pi]), {m, -Infinity, Infinity}];



Generate data:

data = RandomVariate[dis, 10^4];


Match:

FindDistribution[data]


NormalDistribution[1.60001, 1.97652]

Check:

Log[5] // N


1.6094379124341003

So you did not get symbolic Log[5] but at least you got distribution right. This will work in much more complicated cases.

• Possible continuous distributions for TargetFunctions are: BetaDistribution, CauchyDistribution, ChiDistribution, ChiSquareDistribution, ExponentialDistribution, ExtremeValueDistribution, FrechetDistribution, GammaDistribution, GumbelDistribution, HalfNormalDistribution, InverseGaussianDistribution, LaplaceDistribution, LevyDistribution, LogisticDistribution, LogNormalDistribution, MaxwellDistribution, NormalDistribution, ParetoDistribution, RayleighDistribution, StudentTDistribution, UniformDistribution, WeibullDistribution.

• Possible discrete distributions for TargetFunctions are: BenfordDistribution, BinomialDistribution, BorelTannerDistribution, DiscreteUniformDistribution, GeometricDistribution, LogSeriesDistribution, NegativeBinomialDistribution, PascalDistribution, PoissonDistribution, WaringYuleDistribution, ZipfDistribution.

f[m_] = 1/(2*E^((-m + Log[5])^2/8)*Sqrt[2*Pi]);

Integrate[f[m], {m, -Infinity, Infinity}]


1

dist = ProbabilityDistribution[f[m], {m, -Infinity, Infinity}];


Since the integral of f[m] is unity, f[m] does not have to be scaled to be a distribution. A candidate distribution will probably have two parameters and must be defined on the interval {-Infinity, Infinity}. The built-in distributions with two parameters and defined on the interval {-Infinity, Infinity} are

distributions = ToExpression /@ Names["*Distribution"];

twoParamDist = Select[distributions,
FreeQ[PDF[#[a, b], x], #] &&
DistributionDomain[#[a, b]] === Interval[{-Infinity, Infinity}] &] //
Quiet;


The target PDF for comparison is

PDF[dist, m]


Comparing this PDF with the form of those of the built-in two-parameter distributions will greatly limit the candidates

({#, PDF[#[a, b], m]} & /@ twoParamDist) //
Prepend[#, {Distribution, PDF}] & //
Grid[#, Frame -> All] &


The NormalDistribution is the likely distribution

dist2 = NormalDistribution[a, b];

param = Solve[{Mean[dist2] == Mean[dist],
Variance[dist2] == Variance[dist],
DistributionParameterAssumptions[dist2]}, {a, b}]


{{a -> Log[5], b -> 2}}

Verifying that the PDFs are identical

PDF[dist2 /. param[[1]], m] == PDF[dist, m] // Simplify


True

So the target distribution is

dist2 /. param[[1]]


NormalDistribution[Log[5], 2]

• This ought to catch the two-parameter distributions: Select[Names["*Distribution"], (Length["ArgumentsPattern" /. SyntaxInformation[ToExpression[#]]] == 2) &]. After that, you can do further filtering with DistributionDomain[]. Aug 12, 2015 at 3:52
• @J.M. - Your recommendation includes the LandauDistribution. Unfortunately, the PDF for that distribution does not evaluate symbolically. I do not want to have to figure out some kind of workaround for that. Also, similar problem might arise with other distributions when the number of arguments isn't two. Aug 12, 2015 at 5:03
• "Unfortunately, the PDF for that distribution does not evaluate symbolically." - I guess there ought to be a check that the head of the result is not PDF[], but it will admittedly start to look clumsy at that point… Aug 12, 2015 at 5:12
• @J.M.- that is essentially what FreeQ[PDF[#[a, b], x], #]  does. If the PDF doesn't evaluate then the distribution name is still present. Aug 12, 2015 at 5:16
• There are parts of this answer that I like but I think you are using a bit too much human intelligence ;) Yes, once you say that the formula comes from a normal distribution it is not too hard to figure out which normal distribution. But what if the target distribution is not normal and you really don't know what it is. Then what? What similarity metric are you using between the expression for the PDFs and the expression for the target distribution? Aug 12, 2015 at 12:23