Jak obliczyć przedział ufności dla średniej logarytmicznego zestawu danych?

Słyszałem / widziałem w kilku miejscach, że możesz przekształcić zestaw danych w coś o rozkładzie normalnym, biorąc logarytm każdej próbki, oblicz przedział ufności dla przekształconych danych i przekształcaj przedział ufności z powrotem za pomocą operacji odwrotnej (np. podnieś 10 do potęgi dolnej i górnej granicy odpowiednio dla ).$\log_{10}$

Jednak jestem nieco podejrzliwy wobec tej metody, po prostu dlatego, że nie działa ona dla samego środka:$10^{\operatorname{mean}(\log_{10}(X))} \ne \operatorname{mean}(X)$

Jaki jest właściwy sposób to zrobić? Jeśli to nie działa dla samego środka, to jak może działać dla przedziału ufności dla środka?

Istnieje kilka sposobów obliczania przedziałów ufności dla średniej rozkładu logarytmicznego. Przedstawię dwie metody: Bootstrap i prawdopodobieństwo profilu. Przedstawię również dyskusję na temat przeora Jeffreysa.

Bootstrap

Dla MLE

W tym przypadku MLE dla próbki wynosi$(\mu,\sigma)$$(x_1,...,x_n)$

Zatem MLE średniej to . Przez resampling możemy otrzymać próbkę bootstrap z i za pomocą tego, możemy obliczyć kilka bootstrap przedziały ufności. Poniższe kody pokazują, jak je uzyskać.$\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")

Dla próbki średniej

Teraz, biorąc pod uwagę estymator zamiast MLE. Można również rozważyć inny rodzaj estymatorów.$\tilde{\delta}=\bar{x}$

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")

Prawdopodobieństwo profilu

Aby zapoznać się z definicją funkcji wiarygodności i funkcji wiarygodności profilu, patrz . Korzystając z właściwości niezmienniczości prawdopodobieństwa, możemy ponownie sparametryzować w następujący sposób , gdzie a następnie obliczyć numerycznie prawdopodobieństwo profilu .$(\mu,\sigma)\rightarrow(\delta,\sigma)$$\delta=\exp(\mu+\sigma^2/2)$$\delta$

Ta funkcja przyjmuje wartości w ; przedział poziomu ma przybliżoną pewność . Będziemy używać tej właściwości do konstruowania przedziału ufności dla . Poniższe kody pokazują, jak uzyskać ten przedział .$(0,1]$$0.147$ 95 % δ$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$ Bayesian

W tej sekcji przedstawiono alternatywny algorytm oparty na próbkowaniu Metropolis-Hastings i wykorzystaniu wcześniejszych Jeffreysów do obliczania przedziału wiarygodności dla .$\delta$

Przypomnijmy, że w Jeffreys przed o w modelu logarytmiczno-normalnego jest$(\mu,\sigma)$

i że ten przełożony jest niezmienny podczas reparametryzacji. To wcześniejsze jest niewłaściwe, ale tylna część parametrów jest właściwa, jeśli wielkość próbki . Poniższy kod pokazuje, jak uzyskać 95% przedział wiarygodności za pomocą tego modelu Bayesa.$n\geq 2$R

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))

Zauważ, że są bardzo podobne.

Możesz spróbować podejścia bayesowskiego z przeorem Jeffreysa. Powinien dawać przedziały wiarygodności z poprawną właściwością dopasowania częstości: poziom ufności przedziału wiarygodności jest zbliżony do jego poziomu wiarygodności.

 # required package
 library(bayesm)

 # simulated data
 mu <- 0
 sdv <- 1
 y <- exp(rnorm(1000, mean=mu, sd=sdv))

 # model matrix
 X <- model.matrix(log(y)~1)
 # prior parameters
 Theta0 <- c(0)
 A0 <- 0.0001*diag(1)
 nu0 <- 0 # Jeffreys prior for the normal model; set nu0 to 1 for the lognormal model
 sigam0sq <- 0
 # number of simulations
 n.sims <- 5000

 # run posterior simulations
 Data <- list(y=log(y),X=X)
 Prior <- list(betabar=Theta0, A=A0, nu=nu0, ssq=sigam0sq)
 Mcmc <- list(R=n.sims)
 bayesian.reg <- runireg(Data, Prior, Mcmc)
 mu.sims <- t(bayesian.reg$betadraw) # transpose of bayesian.reg$betadraw
 sigmasq.sims <- bayesian.reg$sigmasqdraw

 # posterior simulations of the mean of y: exp(mu+sigma²/2)
 lmean.sims <- exp(mu.sims+sigmasq.sims/2)

 # credibility interval about lmean:
 quantile(lmean.sims, probs = c(0.025, 0.975))

Jednak jestem nieco podejrzliwy wobec tej metody, po prostu dlatego, że nie działa ona dla samego środka: 10mean (log10 (X)) ≠ mean (X)

Masz rację - to wzór na średnią geometryczną, a nie średnią arytmetyczną. Średnia arytmetyczna jest parametrem z rozkładu normalnego i często nie ma większego znaczenia dla danych logarytmicznych. Średnia geometryczna jest odpowiednim parametrem z rozkładu logarytmicznego, jeśli chcesz bardziej sensownie mówić o centralnej tendencji dla twoich danych.

I rzeczywiście obliczalibyśmy CI względem średniej geometrycznej, biorąc logarytmy danych, obliczając średnią i CI jak zwykle i przekształcając wstecz. Masz rację, że tak naprawdę nie chcesz mieszać swoich rozkładów, umieszczając CI dla średniej geometrycznej wokół średniej arytmetycznej ... tak!