Log Normal Distribution

Someone asked me how to generate and compare log-normal distributed data. Here is a short program.

>data=exp(normal(1000)); ...
>d=0.2; {x,y}=histo(data,v=0:d:10);  ...
>plot2d(x,y/1000/d,>bar,fillcolor=lightgray); ...
>plot2d("qnormal(ln(x))/x",epsilon,10,color=red,thickness=3,>add):

test-001

The simulated data can easily be generated by taking the exponential of normal distributed data. The plot above shows how to compute a histogram in EMT. You have to normalize the histogram, so that its components would add to 1, if you went over the complete range of data. Due to the vector v, the data is only counted in between the intervals determined by v. You also have to divide by the lengths of the intervals, since you want to plot the proportion of data in an interval divided by the length of the interval.

The density (pdf) of the log-normal distribution is easy to get due to the substitution t=log(x).

\(\int\limits_{-\infty}^\infty q(t) \, dt = \int\limits_0^\infty q(\log x) \dfrac{dx}{x}\)

The cumulative distribution (cdf) can be computed with a numerical integration as follows.

>function map lognormaldis(x) := integrate("qnormal(ln(x))/x",0,x);
>plot2d("lognormaldis(x)",0,20):

test-002

To compute the inverse, we use the Newton algorithm. This is natural because the derivative of the cumulative distribution is already know.

>function map invlognormaldis(p) ...
$return newton("lognormaldis(x)","qnormal(ln(x))/x",sqrt(p),y=p);
$endfunction
>invlognormaldis(0.9999), lognormaldis(%)
 41.2238299278
 0.9999

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht. Erforderliche Felder sind mit * markiert.

Diese Website verwendet Akismet, um Spam zu reduzieren. Erfahre mehr darüber, wie deine Kommentardaten verarbeitet werden.