# How to use ProbabilityDistribution for a posterior expression [duplicate]

I am performing a Bayesian updating on a Lognormal prior. After combining the prior and likelihood and dividing by the normalization constant I get the following expression as the posterior:

(0.805745 E^(-7.05*10^6 t - 0.255194 (15.576 + Log[t])^2))/t


I now want to use this as a distribution:

pdist = ProbabilityDistribution[posterior, {t, 10^-12, 10^-5}];


So far ok. I can now plot with the distribution

LogLinearPlot[PDF[pdist, t], {t, 10^-11, 10^-6}, PlotRange -> All, FrameLabel -> {"Failure Rate (f/hr)", "PDF"}]


But I cannot take the mean, Mean[pdist], or variance etc, they just return the input. Trying RandomReal[pdist] returns the following:

RandomReal::udist: The specification is not a random distribution recognized by the system


What am I missing in the assignment of the expression to ProbabilityDistribution[].

Clay

• Check out "Create your own distribution" on this page – ssch Oct 17 '13 at 14:04
• @ssch There does not appear to be anything overtly wrong with the OP's set up of the custom distribution. – wolfies Oct 17 '13 at 15:08
• To the OP ... is there a reason your density is numerical/approximate .. rather than setting it up as an exact symbolic algebraic entity? That goes for the bounds of the domain of support too. – wolfies Oct 17 '13 at 15:14

To take a random sample from a distribution you want RandomVariate as opposed to RandomReal, not that it works right away in this case:

posterior = Rationalize[(0.805745 E^(-7.05*10^6 t - 0.255194 (15.576 + Log[t])^2))/t ]
pdist = ProbabilityDistribution[posterior, {t, 10^-12, 10^-5}];
pts = RandomVariate[pdist, 1000];
pts // DeleteDuplicates // Length
(* 1 so it gave the same number a thousand times*)

pts = RandomVariate[pdist, 1000, WorkingPrecision -> 100];
pts // DeleteDuplicates // Length
(* 1000 this seems better *)


Always good to have a look at the quantile plot: But to get Mean, Median, ... to work properly you usually have to define your own distribution, often it is easier to use NIntegrate and FindRoot applied to the definitions of the properties. Together with interpolations of the cdf and inverse cdf this tends to be quite quick.