MC lognormal random distribution is not working as expected

[malcolm-dsider]

The MC lognormal random distribution is not working as expected - an input for HIP-RA such as:

INPUT, Reservoir Porosity, lognormal, 7, 0.5

It should produce a random distribution of numbers centered around a mean of 7, with a standard deviation of 0.5, but instead produces values for Reservoir Porosity in the range of ~300 - 10,000, all of which are invalid for Reservoir Porosity, which has a defined range of 0-100.

---

Note: title/description references to MC/HIP-RA edited for accuracy by @softwareengineerprogrammer

[softwareengineerprogrammer]

Did some digging on this -

1. This appears to have always been the behavior of the lognormal distribution (or it has at least been since the code was brought under source control); see https://github.com/NREL/GEOPHIRES-X/blame/96a1e091fbbd03723270d66a4106936b9e798580/src/geophires_monte_carlo/MC_GeoPHIRES3.py#L195
2. Passing mean/std to np.random.lognormal "draws samples from a log-normal distribution with specified mean, standard deviation, and array shape. Note that the mean and standard deviation are not the values for the distribution itself, but of the underlying normal distribution it is derived from." In other words - the behavior you are seeing is "expected" in the sense that passing 7, 0.5 to lognormal does not produce a distribution whose mean and standard deviation are such.
3. Can you specify the actual desired behavior? If the intention is to produce a lognormal distribution with skew, then the input is insufficiently parameterized and we'd need further information about how you want to specify additional required arguments for skew/moment. (If the desire is to produce a "lognormal" distribution with input mean/std, then this is just equivalent to a normal distribution - see scratch work below)

---

Scratch work showing equivalence between normal and lognormal with mean/std transformed per https://en.wikipedia.org/wiki/Log-normal_distribution#Definitions

    lognormal_vals = []
    mean = 7.0
    std = 0.5
    mu = math.log(mean ** 2 / (np.sqrt(mean ** 2 + std ** 2)))
    sigma = np.sqrt(math.log(1+(std**2)/(mean**2)))

    for i in range(10000):
        lv =np.random.lognormal(mu, sigma)
        lognormal_vals.append(lv)

    def plot(s): # mu, sigma
        #s = np.random.lognormal(mu, sigma, 1000)
        import matplotlib.pyplot as plt
        count, bins, ignored = plt.hist(s, 100, density=True, align='mid')
        x = np.linspace(min(bins), max(bins), 10000)
        pdf = (np.exp(-(np.log(x) - mu) ** 2 / (2 * sigma ** 2))
               / (x * sigma * np.sqrt(2 * np.pi)))
        plt.plot(x, pdf, linewidth=2, color='r')
        plt.axis('tight')
        plt.show()

    plot(lognormal_vals)

[malcolm-dsider]

"Can you specify the actual desired behavior? If the intention is to produce a lognormal distribution with skew, then the input is insufficiently parameterized and we'd need further information about how you want to specify additional required arguments for skew/moment. (If the desire is to produce a "lognormal" distribution with input mean/std, then this is just equivalent to a normal distribution - see scratch work below)"

The expected behavior should be the same as the other random functions - if 7, 0.5 are specified, then it should return a random number close to 7... a number falling within +/- deviation of 0.5 and having a distribution that looks like this:

(see red line)

[softwareengineerprogrammer]

@malcolm-dsider To clarify - the scratch code I included, when run, produces the equivalent plot as your graph, per the equations in the wikipedia page you got it from:

The pdf from the wiki page, however, has distributions with a skew/moment; without skew/moment, the pdf is just equivalent to a normal distribution as you can see. I tagged this issue with question and assigned to you to provide further details on how you want to specify those parameters and how they should behave.