Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

This question already has an answer here:

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[].


share|improve this question

marked as duplicate by ssch, m_goldberg, Sjoerd C. de Vries, Artes, R. M. Oct 18 '13 at 2:06

This question was marked as an exact duplicate of an existing question.

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:

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.

share|improve this answer

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