Генерация логически нормально распределенного случайного числа из среднего, коэффициента вариации

Question

Большинство функций для генерации логически нормально распределенных случайных чисел принимают среднее значение и стандартное отклонение соответствующего нормального распределения в качестве параметров.

Моя проблема в том, что я знаю только среднее значение и коэффициент вариации логнормального распределения. Достаточно просто получить нужные мне параметры для стандартных функций из того, что у меня есть:

Если mu и sigma являются средним и стандартным отклонением соответствующего нормального распределения, мы знаем, что

coeffOfVar^2 = variance / mean^2
             = (exp(sigma^2) - 1) * exp(2*mu + sigma^2) / exp(mu + sigma^2/2)^2
             = exp(sigma^2) - 1

Мы можем изменить это на

sigma = sqrt(log(coeffOfVar^2 + 1))

Мы также знаем, что

mean = exp(mu + sigma^2/2)

Это переставляет на

mu = log(mean) - sigma^2/2

Вот моя реализация R

rlnorm0 <- function(mean, coeffOfVar, n = 1e6)
{
   sigma <- sqrt(log(coeffOfVar^2 + 1))
   mu <- log(mean) - sigma^2 / 2
   rlnorm(n, mu, sigma)
}

Хорошо работает при малых коэффициентах вариации

r1 <- rlnorm0(2, 0.5)
mean(r1)                 # 2.000095
sd(r1) / mean(r1)        # 0.4998437

но не для больших значений

r2 <- rlnorm0(2, 50)
mean(r2)                 # 2.048509
sd(r2) / mean(r2)        # 68.55871

Чтобы проверить, что это не специфичная для R проблема, я переопределил ее в MATLAB. (Использует набор инструментов статистики.)

function y = lognrnd0(mean, coeffOfVar, sizeOut)
if nargin < 3 || isempty(sizeOut)
   sizeOut = [1e6 1];
end
sigma = sqrt(log(coeffOfVar.^2 + 1));
mu = log(mean) - sigma.^2 ./ 2;
y = lognrnd(mu, sigma, sizeOut);
end

r1 = lognrnd0(2, 0.5);
mean(r1)                  % 2.0013
std(r1) ./ mean(r1)       % 0.5008

r2 = lognrnd0(2, 50);
mean(r2)                  % 1.9611
std(r2) ./ mean(r2)       % 22.61

Та же проблема. Вопрос в том, почему это происходит? Это просто, что стандартное отклонение не является устойчивым, когда вариация настолько велика? Или я где-то облажался?

Answer 1

Результаты не удивительны. Для распределений с большим эксцессом ожидаемая дисперсия выборочной дисперсии примерно равна mu4 / N, где mu4 - 4-й момент распределения. Для логнормального значения mu4 экспоненциально зависит от параметра sigma ^ 2, что означает, что для достаточно больших значений sigma ваша выборочная дисперсия будет повсеместно относительно истинной дисперсии. Это именно то, что вы наблюдали. В вашем примере mu4 / N ~ (coeffOfVar ^ 8) / N ~ 50 ^ 8 / 1e6 ~ 4e7.

Для получения ожидаемой дисперсии выборки var. см. http://mathworld.wolfram.com/SampleVarianceDistribution.html. Ниже приведен код, который иллюстрирует идеи более точно. Обратите внимание на большое значение как дисперсии выборочной дисперсии, так и ее теоретического ожидаемого значения, даже для coeffOfVar = 5.

exp.var.of.samp.var <- function(n,mu2,mu4){
  (n-1)*((n-1)*mu4-(n-3)*mu2^2)/n^3
}
mu2.lnorm <- function(mu,sigma){
  (exp(sigma^2)-1)*exp(2*mu+sigma^2)
}
mu4.lnorm <- function(mu,sigma){
  mu2.lnorm(mu,sigma)^2*(exp(4*sigma^2)+2*exp(3*sigma^2)+3*exp(2*sigma^2)-3)
}
exp.var.lnorm.var <- function(n,mu,sigma){
  exp.var.of.samp.var(n,mu2.lnorm(mu,sigma),mu4.lnorm(mu,sigma))
}
exp.var.norm.var <- function(n,mu,sigma){
  exp.var.of.samp.var(n,sigma^2,3*sigma^4)
}    

coeffOfVar <- 5 
mean <- 2
sigma <- sqrt(log(coeffOfVar^2 + 1)) # gives sigma=1.805020
mu <- log(mean) - sigma^2 / 2 # mu=-0.935901
n <- 1e4
m <- 1e4
## Get variance of sample variance for lognormal distribution:
var.trial <- replicate(m,var(rlnorm(n, mu, sigma)))
cat("samp. variance (mean of",m,"trials):",mean(var.trial),"\n")
cat("theor. variance:",mu2.lnorm(mu,sigma),"\n")
cat("variance of the sample var:",var(var.trial),"\n")
cat("expected variance of the sample var:",exp.var.lnorm.var(n,mu,sigma),"\n")
> samp. variance (mean of 10000 trials): 105.7192 
> theor. variance: 100 
> variance of the sample var: 350997.7 
> expected variance of the sample var: 494053.2 
## Do this with normal distribution:
var.trial <- replicate(m,var(rnorm(n, mu, sigma)))
cat("samp. variance (mean of",m,"trials):",mean(var.trial),"\n")
cat("theor. variance:",sigma^2,"\n")
cat("variance of the sample var:",var(var.trial),"\n")
cat("expected variance of the sample var:",exp.var.norm.var(n,mu,sigma),"\n")
> samp. variance (mean of 10000 trials): 3.257944 
> theor. variance: 3.258097 
> variance of the sample var: 0.002166131 
> expected variance of the sample var: 0.002122826