執行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中的平均函數

回答

相關問題