r - 在 R 中实现 1/mean 函数

Question

从正常人群中抽样时，我必须模拟 1/Xbar 的抽样分布。我只想知道我是否正确地开始了我的代码，因为其他一切都取决于此。

MC <- 10000 # Number of samples to simulate

sampling.tau <- function(mu, sigma, sampleSize, MC) {
  tau_hat = c(1:MC)
  for(i in 1:MC)
  {
    mySample <- rnorm(n=sampleSize, mean=mu, sd=sigma)
    tau_hat[i] <- 1/mean(mySample)
  }
}

Answer 1

简短的回答是您走在正确的轨道上。这是一种确认方法。

如果X 是一个分布为 N(μ,σ²) 的随机变量，而 Y 是一个以均值形式形成的随机变量X(例如，X 的 n 个独立样本的总和除以 n)，然后

X ~ N(μ,σ²)

Y ~ N(μ,σ²/n)

您需要来自 Z = 1/Y 分布的样本。一般来说，如果 Y 的密度函数由

给出

Prob(y ≤ Y ≤ y+dy) ≡ f_Y(y), then, if Z = 1/Y

Prob(z ≤ Z ≤ z+dz) ≡ f_Z(z) = (1/z²) × f_Y(1/z)

自

f_Y(y) = √(n/2π) × (1/σ) × exp[-n × (y - μ)²/2σ²]

f_Z(z) = (1/z²) × 1/√2π × (n/√σ) × exp[-n × (1/z - μ)²/2σ²]

所以问题是:您的代码是否会生成分布为 Z 的随机样本？答案可以显示为"is"。

f <- function(z,n,mu=0,sigma=1) 
  (1/z^2)*sqrt(n/(2*pi))*(1/sigma)*exp(-(1/z-mu)^2*(n/(2*sigma^2)))

g <- function(mu, sigma, sampleSize, MC)
  replicate(MC, 1/mean(rnorm(sampleSize, mu, sigma)))

set.seed(1)
hist(g(0,0.1,100,1000),breaks=c(-Inf,seq(-300,300,10),Inf)
     ,xlim=c(-300,300), xlab = "Z",
     main="Histogram of 1/mean(X)", sub="mu=0, sigma=0.1, n=100")
z <- seq(-300,300,1)
lines(z,f(z,100,mu=0,sigma=.1),col="red")