从均值生成对数正态分布的随机数，变异系数
用于生成对数正态分布随机数的大多数函数将相关正态分布的均值和标准差作为参数。
我的问题是我只知道对数正态分布的均值和变异系数。 从我拥有的标准函数中推导出我需要的参数是相当直接的：
如果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
同样的问题。 问题是，为什么会发生这种情况？ 当变化很大时，标准偏差是否不稳定？ 还是我在某个地方搞砸了？
结果并不令人意外。 对于具有大峰度的分布，样本方差的预期方差大致为mu4 / N，其中mu4是分布的第四个时刻。 对于一个对数正态分布，mu4指数地取决于参数sigma ^ 2，这意味着对于足够大的sigma值，您的样本方差将相对于真实方差遍布整个位置。 这正是你所观察到的。 在你的例子中，mu4 / N〜（coeffOfVar ^ 8）/ N〜50 ^ 8 / 1e6〜4e7。
为了推导样本变量的预期方差。 请参阅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
上一篇: Generating lognormally distributed random number from mean, coeff of variation
下一篇: Random sample from given bivariate discrete distribution