Inverse Gaussian with small mean has unreliable sample mean

Question

Consider an inverse Gaussian with a small mean $\mu$ (such as $0.001$) and such that we fix its variance to a larger value $\sigma^2$, say $0.5$. Then if I sample $N$ times from this distribution, I expect the sample mean to be approximately $\mu$, but this does not happen except for wildly large values of $N$.

import scipy as sp

N = 1_000  # we need to use at least 10_000_000, otherwise it is incredibly small, say 1e-6 
mean = 0.001
var = 0.5
lmbda = mean**3 / var

print(f"Sample mean: {sp.stats.invgauss.rvs(mu=(mean/lmbda), loc=0, scale=lmbda, size=N).mean()}")

If you run this script for N equal 1000, 10_000, 1_000_000 you will get around 1e-6 which is in the order of magnitude of the square of the mean.

Why do we get this behavior and how can I avoid this? I.e. can I sample from an inverse gaussian with a small mean and a relatively larger variance without incurring in numerical issues?

I think the reason is due to its incredibly high skewness. If we fix the variance to $\sigma^2$ then we have $\lambda = \mu^3 / \sigma^2$ and using the formula for the skewness of an Inverse Gaussian we get $3\left(\frac{\mu}{\lambda}\right)^{1/2} = 3 \frac{\sigma}{\mu}$. Which means that we get moderate skewness with $\sigma\sim \mu$. In our case, however we have $3\times 0.25 \times 10^{3} = 750$ which is an absurdly large value for skewness — Euler_Salter
– Euler_Salter, Commented Sep 22, 2024 at 21:31
What does "approximately" mean in this context? $\pm 0.1$? $0.01$? Note that the standard deviation is roughly 700x the mean, so the sample standard deviation of the sample mean won't equal the mean until $N \approx 500,000$. — jbowman
– jbowman, Commented Sep 22, 2024 at 21:40
@Euler_Salter - it does not matter much, but the standard deviation is $\sqrt{0.5}$ rather than $0.5^2$ — Henry
– Henry, Commented Sep 22, 2024 at 21:41
The sampling distribution of the mean of iid IG($\mu$,$\lambda$) rvs with mean $\mu$ and variance $\mu^3/\lambda$ has itself an IG($\mu$,$n\lambda$) distribution with mean $\mu$ and variance $\mu^3/n\lambda$ (see here). So for example with $n=1000$, the sample mean has an IG($0.001$, $1000\times (0.001^3/0.5)$) distribution. This distribution is very skewed and the probability that the sample mean is below, say, $0.001/5$ (one fifth of the population mean) is still 0.922. — COOLSerdash
– COOLSerdash, Commented Sep 22, 2024 at 21:42
1. Please consider showing some sample output from your script. 2. Note that the variance of the sample variance (for iid observations) involves the fourth moment of the distribution (and the skewness of that sample variance on the sixth moment). If I did it right, a lower bound on the standard error of the sample variance in very large samples will be on the order of $\sigma^2 |\gamma_1|/\sqrt{n}$ where $\gamma_1$ is the skewness. So with mean $0.001$ and variance $\frac12$, the standard error of the sample variance is $> 1000/\sqrt{n}$ (perhaps much larger, I didn't work its value out). — Glen_b
– Glen_b, Commented Sep 22, 2024 at 22:58

Henry · Accepted Answer · 2024-09-23 10:07:45Z

In your example, the variance of $0.5$ implies a standard deviation of about $0.707$ which is $707$ times the mean of a non-negative random variable; that is a huge coefficient of variation and suggests you are going to need a large sample size to bring the standard error of the mean down to the same order of magnitude as the mean itself.

Things may be slightly worse than that: the skewness is about $2121$, again a very large number, resulting in most of the observations being very small but increasing the probability of seeing extremely large individual observations, so seeing relatively extreme sample means even with moderately large samples sizes.

As an illustration, the first chart below shows the probability of the sample mean being below $\frac12 \mu$ (blue), below $\mu$ (black), and below $\frac32 \mu$ (green) for different sample sizes on the $x$-axis. So it is only for large samples that you are likely to see a sample mean even half of the expectation of the distribution.

produced with the R code

library(statmod)
m <- 0.001
v <- 0.5 
curve(pinvgauss(m*0.5, mean=m, shape=m^3/(v/x)), 
                from=1, to=10^9, log="x", col="blue")  
curve(pinvgauss(m*1.0, mean=m, shape=m^3/(v/x)), 
                from=1, to=10^9, log="x", col="black", add=TRUE)  
curve(pinvgauss(m*1.5, mean=m, shape=m^3/(v/x)), 
                from=1, to=10^9, log="x", col="green", add=TRUE)

while the density curves for the sample means look like the next chart, with $n=10^5$ (red, coefficient of variation $2.24$, skewness $6.71$) and $n=10^6$ (blue, coefficient of variation $0.707$, skewness $2.12$) and $n=10^7$ (green, coefficient of variation $0.224$, skewness $0.671$). As you have observed, it is only with very large sample sizes that something shaped like a normal distribution starts to appear.

produced with the R code

curve(dinvgauss(x, mean=m, shape=m^3/(v/10^5)), 
                from=0, to=m*3, col="red", n=50001)
curve(dinvgauss(x, mean=m, shape=m^3/(v/10^6)), 
                from=0, to=m*3, col="blue", n=50001, add=TRUE)
curve(dinvgauss(x, mean=m, shape=m^3/(v/10^7)), 
                from=0, to=m*3, col="green", n=50001, add=TRUE)
abline(v=m)

For clarity, samples in extremely large individual samples could be replaced by observations. — Richard Hardy
– Richard Hardy, Commented Sep 23, 2024 at 7:04
@RichardHardy Thank you. I have made that change and added a couple of charts — Henry
– Henry, Commented Sep 23, 2024 at 10:08

JimB · Accepted Answer · 2024-09-23 15:35:22Z

The behavior of the sample mean (which is also the maximum likelihood estimator of the mean) is not "unreliable". The behavior is "expected" given the parameters you've chosen (as mentioned by several others). If those are the parameters of interest, you can't avoid that behavior - unless you can show that the issue is a numerical precision issue with the software you're using.

First, looking at the cdf of the sample mean with a sample size of 1,000 is informative. As @CoolSerdash mentions the sample mean also has an inverse gaussian distribution. Using a symbolic algebra system can help examine that issue. Here I'm using Mathematica but I'm sure Maple, MATLAB, and others can do the same.

mu = 1/1000;
variance = 1/2;
lambda = mu^3/variance;
Plot[CDF[InverseGaussianDistribution[mu, 1000 lambda], x], {x, 0, mu},
 PlotRange -> {Automatic, {0, 1}}, Frame -> True, 
 FrameLabel -> {"Sample mean\nwith sample size = 1,000", "Probability density"}]

We see that the probability of a sample mean being less that the mean of the distribution is close to 1:

probMu = CDF[InverseGaussianDistribution[mu, 1000 lambda], mu]

$$\frac{1}{2} \sqrt[250]{e} \left(1-\text{erf}\left(\frac{1}{5 \sqrt{10}}\right)\right)+\frac{1}{2}$$

or

N[CDF[InverseGaussianDistribution[mu, 1000 lambda], mu], 50]
(* 0.96622621731001343285562130481903284097117599744242 *)

You observed that the sample means that you generated are around 1e-6. The probability of a single sample mean being less than 1e-6 is given by

N[CDF[InverseGaussianDistribution[mu, 1000 lambda], 10^-6], 50]
(* 0.15761401956237009972615856026361330363588326211762 *)

The median of the distribution with 1,000 samples is approximately $4.376 \times 10^{-6}$ (here, too, much, much less than the mean).

One could also investigate how the function you're using generates random samples. That could be investigated using arbitrary-precision arithmetic.

Note: While my answer probably sounds like an extended comment, I hope it qualifies as an answer in that it is a response and evidence against the claim that the observed behavior is unreliable.

Stack Exchange Network

Inverse Gaussian with small mean has unreliable sample mean

2 Answers 2

Your Answer

Hot Network Questions

Inverse Gaussian with small mean has unreliable sample mean

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Related

Hot Network Questions