1
$\begingroup$

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

Improve this question
$\endgroup$
3
  • $\begingroup$ Check out "Create your own distribution" on this page $\endgroup$
    – ssch
    Commented Oct 17, 2013 at 14:04
  • $\begingroup$ @ssch There does not appear to be anything overtly wrong with the OP's set up of the custom distribution. $\endgroup$
    – wolfies
    Commented Oct 17, 2013 at 15:08
  • $\begingroup$ 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. $\endgroup$
    – wolfies
    Commented Oct 17, 2013 at 15:14

1 Answer 1

Reset to default
1
$\begingroup$

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:

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.

Improve this answer
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.