13
$\begingroup$

I've fitted a lognormal model using R with a set of data. The resulting parameters were:

meanlog = 4.2991610 
sdlog = 0.5511349

I'd like to transfer this model to Scipy, which I've never used before. Using Scipy, I was able to get a shape and scale of 1 and 3.1626716539637488e+90 -- very different numbers. I've also tried to use the exp of the meanlog and sdlog but continue to get bizarre graph.

I've read every doc I can on scipy and am still confused about what the shape and scale parameters mean in this instance. Would it just make sense to code the function myself? That seems prone to errors though, as I'm new to scipy.

SCIPY Lognormal (BLUE) vs. R Lognormal (RED): Scipy Lognormal (BLUE) vs. R Lognormal (RED)

Any thoughts on what direction to take? The data is fit very well with the R model, by the way, so if it looks like something else in Python, feel free to share.

Thank you!

Update:

I'm running Scipy 0.11

Here's a subset of the data. The actual sample is 38k+, with a mean of 81.53627:

Subset:

x
[60, 170, 137, 138, 81, 140, 78, 46, 1, 168, 138, 148, 145, 35, 82, 126, 66, 147, 88, 106, 80, 54, 83, 13, 102, 54, 134, 34]
numpy.mean(x)
99.071428571428569

Alternatively:

I am working on a function to capture the pdf:

def lognoral(x, mu, sigma):
    a = 1 / (x * sigma * numpy.sqrt(2 * numpy.pi) )
    b = - (numpy.log(x) - mu) ^ 2 / (2 * sigma ^ 2)
    p = a * numpy.exp(b)
    return p

However, this give me numbers the following (I tried several in case I was getting the meaning of sdlog and meanlog mixed up):

>>> lognormal(54,4.2991610, 0.5511349)
0.6994656085799437
 >>> lognormal(54,numpy.exp(4.2991610), 0.5511349)
0.9846125119455129
>>> lognormal(54,numpy.exp(4.2991610), numpy.exp(0.5511349))
0.9302407837304372

Any thoughts?

Update:

rerunning with "UPQuark's" suggestion:

shape, loc, scale (1.0, 50.03445923295007, 19.074457156766517)

The shape of the graph is very similar, however, with the peak happening around 21.

$\endgroup$
4
  • $\begingroup$ This question and answer may help: stackoverflow.com/questions/8747761/… $\endgroup$
    – jbowman
    Jul 25, 2012 at 21:07
  • $\begingroup$ Thank you, I found that and learned the "fitting" with the lognormal. However, my questions is why would I get such different distributions? $\endgroup$ Jul 25, 2012 at 21:15
  • $\begingroup$ Are you using SciPy 0.9? Also, could you post your data, or a subset thereof? $\endgroup$
    – jbowman
    Jul 25, 2012 at 21:18
  • $\begingroup$ Updated! It's Scipy 0.11 by the way. So the bugs I've read about shouldn't be relevant ;) $\endgroup$ Jul 25, 2012 at 21:37

4 Answers 4

16
$\begingroup$

I fought my way through the source code, to arrive at the following interpretation of the scipy lognormal routine.

$\frac{x-\text{loc}}{\text{scale}} \sim \text{Lognormal}(\sigma)$

where $\sigma$ is the "shape" parameter.

The equivalence between scipy parameters and R parameter is as follows:

loc - No equivalent, this gets subtracted from your data so that 0 becomes the infimum of the range of the data.

scale - $\exp{\mu}$, where $\mu$ is the mean of the log of the variate. (When fitting, typically you'd use the sample mean of the log of the data.)

shape - the standard deviation of the log of the variate.

I called lognorm.pdf(x, 0.55, 0, numpy.exp(4.29)) where the arguments are (x, shape, loc, scale) respectively, and generated the following values:

x      pdf
10     0.000106
20     0.002275
30     0.006552
40     0.009979
50     0.114557
60     0.113479
70     0.103327
80     0.008941
90     0.007494
100    0.006155

which seem to match pretty well with your R curve.

$\endgroup$
1
  • $\begingroup$ Thank you, @JBowman, that is exactly the explanation I needed and the output is precisely my distribution. $\endgroup$ Jul 26, 2012 at 17:13
8
$\begingroup$

The lognormal distribution in SciPy fits in to the general framework for all distributions in SciPy. They all have a scale and location keyword (which default to 0 and 1 if not explicitly provided). This allows all distributions to be shifted and scaled from their normalized specification with clear implications to the statistics of the distribution. The distributions typically have one or more "shape" parameters as well (though some, like the normal distribution, do not need any additional parameters).

While this general approach nicely unifies all the distributions, for lognormal it can create some confusion because of the way other packages define the parameters. Still, it is very simple to match any lognormal distribution if you meanlog (the mean of the underlying distribution) and sdlog (the standard deviation of the underlying distribution).

First, make sure you set the location parameter to 0. Then, set the shape parameter to the value of sdlog. Finally, the set the scale parameter to math.exp(meanlog). Thus, rv = scipy.stats.lognorm(0.5511349, scale=math.exp(4.2991610)) will create a distribution object whose pdf matches your R-generated curve exactly. As x = numpy.linspace(0,180,1000); plot(x, rv.pdf(x)) will verify.

Basically, the SciPy lognormal distribution is a generalization of the standard lognormal distribution which matches the standard exactly when setting the location parameter to 0.

When fitting data with the .fit method, you can also use keywords, f0..fn, floc, and fshape to hold fixed any of the shape, location, and/or scale parameters and only fit over the other variables. For the lognormal distribution this is very useful as usually you know the location parameter should be fixed to 0. Thus, scipy.stats.lognorm.fit(dataset, floc=0) will always return the location parameter as 0 and only vary the other shape and scale parameters.

$\endgroup$
3
$\begingroup$

Scipy lognormal fit returns shape, location, and scale. I just ran the following on an array of sample price data:

shape, loc, scale = st.lognorm.fit(d_in["price"])

This gives me reasonable estimates 1.0, 0.09, 0.86, and when you plot it, you should take into account all three parameters.

The shape parameter is the standard deviation of the underlying normal distribution, and the scale is the exponential of the mean of the normal.

Hope this helps.

$\endgroup$
5
  • $\begingroup$ Thank you for replying! Once I have these values (loc, scale, shape), I try to find the pdf(x) for each x that I care about (here it is values from 0 to 180, exclusive). scipy.stats.lognorm.pdf(i, loc, scale, shape) However, plotting these I get the plot above. $\endgroup$ Jul 25, 2012 at 21:40
  • $\begingroup$ OK I saw you mentioning shape and scale only, that's why I mentioned there are three parameters returned by default from fit(). You also said you're confused about what the shape and scale parameters mean, and I tried to address that. I've never had the lognormal fit return absurd values like in your case though, what is the location parameter? $\endgroup$
    – upquark
    Jul 25, 2012 at 21:54
  • $\begingroup$ Just updated the question to answer that. Thanks for thinking about this. $\endgroup$ Jul 25, 2012 at 22:00
  • $\begingroup$ Call scipy.stats.lognorm.pdf(x, shape, loc, scale) instead of scipy.stats.lognorm.pdf(i, loc, scale, shape). $\endgroup$
    – upquark
    Jul 25, 2012 at 22:10
  • $\begingroup$ Thanks, upquark, I've done that as well with similar results. The whole shape of the graph continues to be very different from the expected results giving in R. Looks like a totally different distribution than the one in R, actually. $\endgroup$ Jul 26, 2012 at 13:22
1
$\begingroup$

Seems like the distribution in Scipy for the lognormal is not the same as in R, or generally, not the same as the distribution I am familiar with. John D Cook has touched on this: http://www.johndcook.com/blog/2010/02/03/statistical-distributions-in-scipy/ http://www.johndcook.com/distributions_scipy.html

However, I haven't found anything conclusive on how to use a lognormal density function in Python. If anyone would like to add to this, please feel free.

My solution so far is to use the lognormal pdf evaluated at 0 to 180 (exclusive), and used as a dictionary in the python script.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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