Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

Clay

share|improve this 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
add comment

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

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

1 Answer

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
add comment

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