如何計算對數正態數據集均值的置信區間？
我在幾個地方聽說/看到您可以通過取每個樣本的對數將數據集轉換為正態分佈的數據，計算轉換數據的置信區間，然後使用逆運算將置信區間轉換回來（例如，分別提高 10 的下限和上限的冪，對於）。
但是，我對這種方法有點懷疑，僅僅是因為它不適用於均值本身：
這樣做的正確方法是什麼？如果它不適用於平均值本身，它怎麼可能適用於平均值的置信區間？
有幾種方法可以計算對數正態分佈均值的置信區間。我將介紹兩種方法：Bootstrap 和 Profile 可能性。我還將介紹關於 Jeffreys 之前的討論。
引導程序
對於 MLE
在這種情況下，MLE $ (\mu,\sigma) $ 樣品 $ (x_1,…,x_n) $ 是
$$ \hat\mu= \dfrac{1}{n}\sum_{j=1}^n\log(x_j);,,,\hat\sigma^2=\dfrac{1}{n}\sum_{j=1}^n(\log(x_j)-\hat\mu)^2. $$
那麼，均值的 MLE 為 $ \hat\delta=\exp(\hat\mu+\hat\sigma^2/2) $ . 通過重採樣，我們可以獲得一個自舉樣本 $ \hat\delta $ 並且，使用它，我們可以計算幾個引導置信區間。以下
R代碼顯示瞭如何獲取這些。rm(list=ls()) library(boot) set.seed(1) # Simulated data data0 = exp(rnorm(100)) # Statistic (MLE) mle = function(dat){ m = mean(log(dat)) s = mean((log(dat)-m)^2) return(exp(m+s/2)) } # Bootstrap boots.out = boot(data=data0, statistic=function(d, ind){mle(d[ind])}, R = 10000) plot(density(boots.out$t)) # 4 types of Bootstrap confidence intervals boot.ci(boots.out, conf = 0.95, type = "all")對於樣本均值
現在，考慮估計器 $ \tilde{\delta}=\bar{x} $ 而不是 MLE。也可以考慮其他類型的估計器。
rm(list=ls()) library(boot) set.seed(1) # Simulated data data0 = exp(rnorm(100)) # Statistic (MLE) samp.mean = function(dat) return(mean(dat)) # Bootstrap boots.out = boot(data=data0, statistic=function(d, ind){samp.mean(d[ind])}, R = 10000) plot(density(boots.out$t)) # 4 types of Bootstrap confidence intervals boot.ci(boots.out, conf = 0.95, type = "all")輪廓可能性
有關似然函數和輪廓似然函數的定義，請參閱。使用似然性的不變性，我們可以重新參數化如下 $ (\mu,\sigma)\rightarrow(\delta,\sigma) $ ， 在哪裡 $ \delta=\exp(\mu+\sigma^2/2) $ 然後數值計算輪廓似然 $ \delta $ .
$$ R_p(\delta)=\dfrac{\sup_{\sigma}{\mathcal L}(\delta,\sigma)}{\sup_{\delta,\sigma}{\mathcal L}(\delta,\sigma)}. $$
這個函數取值 $ (0,1] $ ; 水平區間 $ 0.147 $ 有一個近似的置信度 $ 95% $ . 我們將使用這個屬性來構建一個置信區間 $ \delta $ . 以下
R代碼顯示瞭如何獲取此間隔。set.seed(1) # Simulated data data0 = exp(rnorm(100)) # Log likelihood ll = function(mu,sigma) return( sum(log(dlnorm(data0,mu,sigma)))) # Profile likelihood Rp = function(delta){ temp = function(sigma) return( sum(log(dlnorm(data0,log(delta)-0.5*sigma^2,sigma)) )) max=exp(optimize(temp,c(0.25,1.5),maximum=TRUE)$objective -ll(mean(log(data0)),sqrt(mean((log(data0)-mean(log(data0)))^2)))) return(max) } vec = seq(1.2,2.5,0.001) rvec = lapply(vec,Rp) plot(vec,rvec,type="l") # Profile confidence intervals tr = function(delta) return(Rp(delta)-0.147) c(uniroot(tr,c(1.2,1.6))$root,uniroot(tr,c(2,2.3))$root)$ \star $ 貝葉斯
在本節中，一種基於 Metropolis-Hastings 抽樣和使用 Jeffreys 先驗的替代算法，用於計算可信區間 $ \delta $ 被呈現。
回想一下Jeffreys 之前的 $ (\mu,\sigma) $ 在對數正態模型中是
$$ \pi(\mu,\sigma)\propto \sigma^{-2}, $$
並且這個先驗在重新參數化下是不變的。這個先驗是不合適的，但是如果樣本量大，參數的後驗是合適的 $ n\geq 2 $ . 以下
R代碼顯示瞭如何使用此貝葉斯模型獲得 95% 的可信區間。library(mcmc) set.seed(1) # Simulated data data0 = exp(rnorm(100)) # Log posterior lp = function(par){ if(par[2]>0) return( sum(log(dlnorm(data0,par[1],par[2]))) - 2*log(par[2])) else return(-Inf) } # Metropolis-Hastings NMH = 260000 out = metrop(lp, scale = 0.175, initial = c(0.1,0.8), nbatch = NMH) #Acceptance rate out$acc deltap = exp( out$batch[,1][seq(10000,NMH,25)] + 0.5*(out$batch[,2][seq(10000,NMH,25)])^2 ) plot(density(deltap)) # 95% credibility interval c(quantile(deltap,0.025),quantile(deltap,0.975))請注意，它們非常相似。