执行1/R中的平均函数

我必须模拟1/Xbar从正常人群采样时的采样分布。我只是想知道是否我的代码是正确的，因为一切都依赖于此。执行1/R中的平均函数

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

来源

2014-02-11 user3295513

首先注意到，有一个在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

来源

2014-02-11 06:15:49 tonytonov

简短的回答是，你在正确的轨道上。这是一种确认它的方法。

如果X是分布式如N的随机变量（μ，σ ），并Y是形成为X平均的随机变量（例如，X的n个独立样本的和，除以n），然后

X〜N（μ，σ ）

Y〜N（μ，σ/n）

你想从Z = 1/Y的分布样本。通常，如果为Y密度函数由

习题（Y ≤ý≤ Y + DY）≡˚F_Ý（y）的给定的，那么，如果Z = 1/Y

习题（Z ≤ž≤ Z + DZ）≡˚F_ž（Z）=（1/Z 2 ）×˚F_Ý（1/Z）

由于

˚F_Ý（Y）= √（N/2 π）×（1/σ）× EXP [-n ×（Y - μ）/2 σ ]

˚F_ž（Z）=（1/Z 2 ） × 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")

来源

2014-02-11 15:57:49 jlhoward

非常漂亮的确。 –

如果你想的1/mean(X)抽样分布，其中X是正常您可以通过认识到如果X已经意味着mu节省自己的时间很多，SD sigma，然后抽样分布N的平均值为正常，平均值为mu和sd sigma/sqrt(N)，所以：

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

应该多效率，并提供可比较的结果（很显然你应该仔细检查它针对蛮力解决方案，为自己...）

来源

2014-02-11 16:28:06

执行1/R中的平均函数

回答

相关问题