0

Here is a probability density function of a lognormal distribution:

from scipy.stats import lognorm
def f(x): return lognorm.pdf(x, s=0.2, loc=0, scale=np.exp(10))

This function has very small y values (max ~ 1E-5) and distributes over x value ~1E5. We know that the integral of a PDF should be 1, but when using the following codes to directly calculate integral, the answer is round 1E-66 since the computation accuracy is not enough.

from scipy.integrate import quad
import pandas as pd
ans, err = quad(f, -np.inf, np.inf)

Could you kindly help me to correctly calculate an integral like this? Thank you.

Yiyin
  • 75
  • 2
  • 9

1 Answers1

2

The values that you are using correspond to the underlying normal distribution having mean mu = 10 and standard deviation sigma = 0.2. With those values, the mode of the distribution (i.e. the location of the maximum of the PDF) is at exp(mu - sigma**2) = 21162.795717500194. The function quad works pretty well, but it can be fooled. In this case, apparently quad only samples the function where the values are extremely small--it never "sees" the higher values way out around 20000.

You can fix this by computing the integral over two intervals, say [0, mode] and [mode, np.inf]. (There is no need to compute the integral over the negative axis, since the PDF is 0 there.)

For example, this script prints 1.0000000000000004

import numpy as np
from scipy.stats import lognorm
from scipy.integrate import quad


def f(x, mu=0, sigma=1):
    return lognorm.pdf(x, s=sigma, loc=0, scale=np.exp(mu))


mu = 10
sigma = 0.2

mode = np.exp(mu - sigma**2)

ans1, err1 = quad(f, 0, mode, args=(mu, sigma))
ans2, err2 = quad(f, mode, np.inf, args=(mu, sigma))

integral = ans1 + ans2
print(integral)
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • Thank you very much for this helpful answer! May I also learn generally to which value will `quad` sample? How is this determined or at what kind of case I have to split the integral to two intervals? Thanks again. – Yiyin Sep 08 '20 at 08:57