在 R 中实现 1/mean 函数答案

【问题标题】：Implementing 1/mean function in R在 R 中实现 1/mean 函数
【发布时间】：2023-03-16 18:29:01
【问题描述】：

从正常人群中抽样时，我必须模拟 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)
  }
}

【问题讨论】：

标签： r function statistics

【解决方案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)，那么，如果 Z = 1/Y

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

自从

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

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

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

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

【讨论】：

【解决方案2】：

首先注意，c(1:MC) 中不需要，只使用1:MC，这已经是一个向量了。第二个注意事项，您不返回tau_hat。第三点，声明tau_hat <- numeric(MC) 可能是更好的方法：无论如何你都在覆盖它。

除此之外，一切看起来都不错。我会稍微修改你的代码以避免循环：

sampling.tau.2 <- function(mu, sigma, sampleSize, MC) {
  replicate(MC, 1/mean(rnorm(sampleSize, mu, sigma)))
}

sampling.tau.2(10, 1, 100, 5)
# values should be close to 1/mu = 1/10
[1] 0.09808410 0.10000718 0.09870573 0.09952546 0.09843164

【讨论】：

【解决方案3】：

如果您想要1/mean(X) 的抽样分布，其中X 是正常的，您可以通过认识到如果X 具有平均mu，sd @987654325，可以节省大量时间 @，则N样本均值的采样分布为Normal，均值为mu，sd为sigma/sqrt(N)，所以：

sampling.tau.3 <- function(mu, sigma, sampleSize, MC) {
   1/rnorm(MC, mu, sigma/sqrt(sampleSize))
}

应该快得多并给出可比较的结果（您显然应该根据自己的蛮力解决方案仔细检查它......）

【讨论】：