6
$\begingroup$

I'd like to generate realistic sample data in javascript for monetary donations. I've chosen to model them as following the log normal distribution.

I think a fairly intuitive way for a lay person to specify the parameters of a function that generates this data is by giving a median and 99th percentile. To put it into every day words, "the average (median) donation is £50 and 1 in 100 donations are over £1000".

I'm trying to write a javascript function that takes the median and 99th percentile of a log normal distribution and outputs the mu and sigma of that distribution so I can use these as inputs to d3.randomLogNormal(mu, sigma).

I've been able to write this function for the normal distribution as follows:

import {randomNormal} from "d3";


function normalGenerator(p50, p99) {
  
  const mu = p50;
  // The 99th percentile is always 2.3263 standard deviations away from the median.
  const sigma = (p99 - p50) / 2.3263;
  
  return randomNormal(mu, sigma);
}

But I am stuck doing this for the log normal distribution. Here is where I have got to so far:

import {randomLogNormal} from "d3";


function logNormalGenerator(p50, p99) {
  
  const mu = Math.log(p50);
  const sigma = // Some function of p50 and p99
  
  return randomLogNormal(mu, sigma);
}

I am fairly sure that mu should be Math.log(p50) but cannot work out how to express sigma in terms of p50 and p99 (despite having read the Wikipedia log normal article and various other SE answers, and attempting some paper and pencil efforts to 'transpose' between normal and log normal).

I've tried combinations of things like (Math.log(p99)-Math.log(p50)) / 2.3263 but at this point, it feels like I am just punching in random combinations.

Perhaps what I need to do is 're-arrange the CDF of the log normal distribution to get sigma on the left. I can then substitute in mu=p50 and x=p99'. If that is the case, then I think that doing so is beyond my maths skills, though I am wondering if there is a shortcut similar to what I did for the normal distribution. In any case, I am feeling a little bit lost with the maths!

If it is helpful, here is an an observable notebook https://observablehq.com/d/ac3ebdbc2ab4af78 that checks these functions by generating a large amount of samples and checking the p50 and p99 of the generated data.

Would be great to get an answer to this and if possible an explanation of the maths behind it for someone who is a bit rusty with logs and has close to zero experience with the log normal distribution :)

Additional info:

I came across https://distribution-explorer.github.io/continuous/lognormal.html#pdf-and-cdf-plots which has a GUI for doing what I want to do (though not in javascript).

If you put it in 'Quantile setter mode' and plug in the following:

  • lower y: 50
  • lower quantile: .5
  • upper y: 1000
  • upper quantile: .99

It spits out μ = 3.912, σ = 1.288.

And when I plug those numbers into my Observable notebook, I do get samples that fit (I've added a row to the table 'Correct log normal' to illustrate).

So essentially, what I am looking for is a JS function to calculate σ from the above four values.

Final thoughts:

After finding a solution to the original question, I came to the conclusion that hard coding a sigma was the best way to go.

When you allow people to specify p50 and p99 based on their intuition, they often supply arguments that generate unrealistic sigmas / shapes, e.g. an unrealistically high sigma that results in large number of values approaching 0. Hard coding sigma as 1 generates some fairly reasonable values.

I also decided to allow them to specify the mean rather than the median since mean is more intuitive to most.

This led to the following function:

function logNormalGenerator(mean) {
  const mu = Math.log(mean) - 1.6487;
  return d3.randomLogNormal(mu, 1);
}

Or if you don't want to hard code sigma (and keep it reasonable via some other method):

function logNormalGenerator(mean, sigma) {
  const mu = Math.log(mean) - ((sigma**2)/2);
  return d3.randomLogNormal(mu, sigma);
}
$\endgroup$
9
  • 1
    $\begingroup$ The log normal distribution has two parameters identifying it, why can the users not specify these directly which would avoid all this? I would imagine that specifying location and scale is easier than specifying the median and 99-th percentile. $\endgroup$ Commented Jun 17 at 10:03
  • 1
    $\begingroup$ Good question :) A lay person would not have an idea of what mu and sigma would translate to in terms of the size of donations. e.g. I wouldn't expect them to know that mu=5 and sigma=1.5 means average donation of £145.45 and 1 in 100 over £5,365.411. Where as they would intuitively know what average donation of £50 and 1 in 100 donations are over £1000" would look like. $\endgroup$ Commented Jun 17 at 10:50
  • 1
    $\begingroup$ The problem with specifying only "one side" of the distribution is that the lognormal distribution is asymmetric and so you have no idea what you will get on the "other side", this only works for symmetric distributions. And if you specify both sides you might get nonsense results as these might be far off from any theoretical distribution. $\endgroup$ Commented Jun 17 at 10:54
  • 1
    $\begingroup$ If you have a function for the normal, just pass the logs of the lognormal quantiles to it. $\endgroup$
    – Glen_b
    Commented Jun 17 at 16:34
  • 1
    $\begingroup$ @michaelmcandrew It should be easy enough: you have $\mu+\frac12 \sigma^2 = \log(lognormalmean)$ and $\mu+2.3263 \sigma = \log(lognormalp99)$ so can substitute to get a quadratic equation for $\mu$ and thus find $\sigma$. $\endgroup$
    – Henry
    Commented Jun 19 at 8:58

2 Answers 2

4
$\begingroup$

If $X \sim \operatorname{Normal}(\mu, \sigma^2)$ then $Y=e^X$ has a log-normal distribution with the same parameters and $X=\log(Y)$. So you simply need to take logarithms and exponentiation appropriately.

You presumably want something like:

function logNormalGenerator(p50, p99) {
  
  const mu = Math.log(p50);
  const sigma = (Math.log(p99) - Math.log(p50)) / 2.3263;
  
  return randomLogNormal(mu, sigma);
}

to match your earlier normalGenerator function.

Alternatively, you could use that function directly with logarithms and exponentiation as in

function logNormalGenerator(p50, p99) {
  
  return Math.exp(normalGenerator(Math.log(p50), Math.log(p99)));
}
$\endgroup$
1
  • $\begingroup$ you are right - I was missing some brackets around the top half of the sigma calculation in the function in the notebook. Adding them in as per your function gives the answers I was looking for. $\endgroup$ Commented Jun 18 at 11:01
6
$\begingroup$

A log-normal variable is normally distributed in the log scale, so the identity holds in the log scale.

To set the scene, let's review the normal case first. The $p^{th}$ quantile $Q(p)$ is calculated like so:

$$ Q(p)=\mu+\sigma\sqrt 2\ \text{erf}^{-1}(2p-1) $$

You can trivially verify that $Q(0.5)=\mu$ since then $\text{erf}^{-1}(0)=0$, so your current calculation is $Q(0.99)-Q(0.5)=\mu+\sigma\sqrt 2 \text{erf}^{-1}(0.98)-\mu=\sigma\sqrt 2\ \text{erf}^{-1}(0.98)$.

The log-normal quantile is very similar, except that it's now the log scale that follows the above formula:

$$ Q(p)=\text{exp}\big(\mu+\sigma\sqrt 2\ \text{erf}^{-1}(2p-1)\big) $$

As before, $Q(0.5)=\text{exp}(\mu)$. Because of the exponentiation, which is inherently a multiplicative operation, it's much easier to express the original scale as a ratio rather than a difference:

$$ \frac{Q(0.99)}{Q(0.5)}=\frac{\text{exp}\big(\mu+\sigma\sqrt 2\ \text{erf}^{-1}(0.98)\big)}{\text{exp}\big(\mu\big)}=\text{exp}\big(\sigma\sqrt 2\ \text{erf}^{-1}(0.98)\big) $$

Simplifying leads to:

$$ \frac{\log\big(Q(0.99)/Q(0.5)\big)}{\sqrt 2\ \text{erf}^{-1}(0.98)}=\sigma $$

An implementation in R could look like so:

## Inverse error function (erf^-1)
erfi <- function(x) qnorm((1+x)/2)/sqrt(2)

## Convert a median and another quantile to normal parameters
q2norm <- function(median, p2, p=.99) {
  sig <- abs(p2-median) / (sqrt(2) * erfi(2*p-1))
  setNames(c(median, sig), c("mu", "sig"))
}

x1 <- rnorm(1E6, mean=1.23, sd=0.56)
q1 <- quantile(x1, probs = c(.5, .95))
q2norm(q1[1], q1[2], p=.95)
#>     mu    sig
#> 1.2306 0.5596


## The same for log-normal
q2lnorm <- function(median, p2, p=.99) {
  mu <- log(median)
  sig <- log(p2/median)/(sqrt(2) * erfi(2*p-1))
  setNames(c(mu, sig), c("mu", "sig"))
}

x2 <- rlnorm(1E7, 3.12, 0.333)
q2 <- quantile(x2, probs = c(.5, .9))
q2lnorm(q2[1], q2[2], p=.9)
#>     mu    sig
#> 3.1197 0.3329

Finally, I do want to stress that extreme sample percentiles from even a theoretically perfect log-normal can be quite unstable, see for example this discussion. There's a reason why I used millions of samples above...

$\endgroup$
5
  • 2
    $\begingroup$ A little jargon heavy for $\exp(X) < q \iff X < \log(q)$, but your final point is really nice +1. Just want to point out that your final formula for $\sigma$ is equivalent to OPs (Math.log(p99)-Math.log(p50)) / 2.3263. Please correct me if I'm wrong. $\endgroup$ Commented Jun 17 at 12:16
  • 2
    $\begingroup$ @LukasLohse your last point is indeed true, I tried to expand a bit on why that is rather than just give the calculation, because otherwise this would be a programming question :o) $\endgroup$
    – PBulls
    Commented Jun 17 at 12:37
  • $\begingroup$ Thanks for the answer @LukasLohse. I am not familiar with R but it looks like the key line is the one that sets sig in the q2lnorm function and that I'll need to find an inverse error function in Javascript to do the same. This one looks hopeful: github.com/stdlib-js/math-base-special-erfinv $\endgroup$ Commented Jun 18 at 8:22
  • $\begingroup$ @michaelmcandrew No you don't! If you stick with 99%-quantile then sqrt(2) * erfi(2*p-1) = 2.3263. $\endgroup$ Commented Jun 18 at 9:00
  • $\begingroup$ OK - thanks for clarifying @LukasLohse! $\endgroup$ Commented Jun 18 at 11:02

Your Answer

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

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