私には関係があります:y = a + b + c
a、b、cの平均と標準偏差があり、モンテカルロシミュレーションによってこれからyの確率分布を取得したいと思います。
これを行うために使用できる関数やパッケージ、または簡単な方法はありますか?
私には関係があります:y = a + b + c
a、b、cの平均と標準偏差があり、モンテカルロシミュレーションによってこれからyの確率分布を取得したいと思います。
これを行うために使用できる関数やパッケージ、または簡単な方法はありますか?
入力a、b、cは平均と標準偏差で定義できると言うので、入力a、b、cが正規分布していると仮定していると思います。その場合は、特別なパッケージがなくても、これをかなり高速に実行できます。
mu.a=33
mu.b=32
mu.c=13
sigma.a=22
sigma.b=22
sigma.c=222
n= a.large.number=10^5
a=rnorm(n,mu.a,sigma.a)
b=rnorm(n,mu.b,sigma.b)
c=rnorm(n,mu.c,sigma.c)
y=a+b+c
plot(density(y))
mean(y)
sd(y)
y
、a
、b
およびについて行っているすべての仮定に注意してくださいc
。y の平均の標本分散を計算するなど、より複雑なことをしたい場合。次に、この手順を何度も実行して平均を収集し、プロットします。
mysimfun=function(n,mu,sigma,stat.you.want='mean')
# mu is length 3 and sigma is too.
{
n= a.large.number=10^5
a=rnorm(n,mu[1],sigma[1])
b=rnorm(n,mu[2],sigma[2])
c=rnorm(n,mu[3],sigma[3])
y=a+b+c
plot(density(y))
return(ifelse(stat.you.want=='mean',mean(y),sd(y))
}
mu=c(mu.a,my.b,mu.c)
sigma=c(sigma.a,sigma.b,sigma.c)
mi=rep(NA,100)
次に、ある種のループで実行します。
for(i in 1:100) {mi[i]=mysimfun(10,mu,sigma,stat.you.want='mean') }
par(mfrow=c(2,1)
hist(mi)
plot(density(mi))
mean(mi)
sd(mi)
2 つのアプローチがあります: モンテカルロが意味するブートストラップ、または経験的分布から推定値を構築するよりも理論に興味がある場合は、「distr」パッケージとその友人である「distrSim」および「distrTEst」です。
require(boot)
ax <- rnorm(100); bx<-runif(100); cx<- rexp(100)
dat <- data.frame(ax=ax,bx=bx,cx=cx)
boot(dat, function(d){ with(d, mean(ax+bx+cx) )}, R=1000, sim="parametric")
boot(dat, function(d){ with(d, sd(ax+bx+cx) )}, R=1000, sim="parametric")