CODEWITHSUNDEEP
pythonscipy
Bobesh
Bobesh

Reputation: 1207

Python lognormal density differs from analytical solution

I have been testing some functions in python (to make sure I understand functions well) and I am confused because I have different results.

I was testing stats.lognorm.pdf form scipy. This function should return the same results with x,shape, scale, loc = 0 in the following code:

val1 = (1/(x*shape*math.sqrt(2*math.pi)))*math.exp(-((math.log(x)-math.log(scale)**2)/(2*math.pow(shape,2))))
val2 = stats.lognorm.pdf(x, shape, 0, scale)    #I expect that val1 == val2

When I try it with some fine numbers it looks fine.

x = 1

scale = 1 #log(shape) = u => u=0 
shape = 0.25

then 

val1 = 1.5957691216057308

val2 = 1.59576912161

but when I set 

shape = 0.8
scale = 25.16
x = 23

the results differ a lot

val1 = 6.33367993244142
val2 = 0.0215455972263

Why is this happening? Is something wrong with my code?

Upvotes: 1

Views: 67

Answers (2)

ymmx
ymmx

Reputation: 4967

I think you misread the documentation

The probability density function for lognorm is:

lognorm.pdf(x, s) = 1 / (s*x*sqrt(2*pi)) * exp(-1/2*(log(x)/s)**2)

for x > 0, s > 0.

lognorm takes s as a shape parameter.

The probability density above is defined in the “standardized” form. To
shift and/or scale the distribution use the loc and scale parameters. 
Specifically, lognorm.pdf(x, s, loc, scale) is identically equivalent to 
lognorm.pdf(y, s) / scale with y = (x - loc) / scale.

If we follow the step we have:

x = 23.
shape = 0.8
scale = 25.16
loc = 0.


xp = (x - loc) / scale
val1 = 1. / (shape*xp*math.sqrt(2.*math.pi)) * math.exp(-1./2.*(math.log(xp)/shape)**2)
val1 = (val1) / scale

print(val1)
val2 = stats.lognorm.pdf(x, shape, 0, scale)    #I expect that val1 == val2
print(val2)

which give : 

0.02154559722626566
0.0215455972263

Upvotes: 1

Jaime
Jaime

Reputation: 67437

Your val1 is wrong, you have the **2 in the exponent inside, rather than outside, the parenthesis. If you try with this:

val1 = (1 / (x * shape * math.sqrt(2 * math.pi))) * math.exp(
      -(math.log(x) - math.log(scale))**2 / (2 * math.pow(shape, 2)))

everything should work as expected.

There may be a lesson to learn here on why PEP8 insists on properly formatting and spacing your code, as it makes spotting errors like this easier.

Upvotes: 1

Related Questions