2016-10-05 12 views
0

私はRには比較的新しいので、R2OpenBUGSを使用してMCMCシミュレーションを行っています。私のモデルはうまく機能し、効率的に複数の反復をしようとしています。これは複数の反復のための私のコードです。レプリケート関数からの値の抽出

mcmc.s=function(n1,m1,k1,al1,be1,R.x,n2,m2,k2,al2,be2,R.y) 
{ 
    tq.x=replicate(n1, aa(k1,al1,be1)) #produces data 
    tq.y=replicate(n2, aa(k2,al2,be2)) #produces data 

    x=pfF(n1,m1,k1,tq.x) #extracts samples from data 
    y=pfF(n2,m2,k2,tq.y) #extracts samples from data 

    x=c(x) #Need data inputed as a list 
    y=c(y) #Need data inputed as a list 

    datatemp<- list("x","y","R.x","k1","m1","R.y","k2","m2") #these are just my inputs 
    aaa<-rgamma(1, shape=4.5,rate=52) #prior 
    bbb<-rgamma(1, shape=4.5,rate=52) #prior 
    ccc<-rgamma(1, shape=20.5,rate=3.5) #prior 

    initstemp<-function(){list(a=aaa,b=bbb,c=ccc)} #assigns priors to parameters in model 
    bugstemp = bugs(data=datatemp,inits=initstemp,parameters.to.save=c("a","b","c"),n.iter= 50000, model.file=mtemp, 
        n.chains=3, n.burnin=40000,n.thin=1,codaPkg =F, debug=F) 
    out<-bugstemp$summary 
    acalc<-out["a","mean"];bcalc<-out["b","mean"];ccalc<-out["c","mean"] #line not necessary 
    return(out) 
} 

iter=2 
out=replicate(iter,mcmc.s(n1,m1,k1,al1,be1,R.x,n2,m2,k2,al2,be2,R.y)) 
out 

結果は次のようになります。

, , 1 

       mean   sd  2.5%  25%  50%  75%  97.5%  Rhat n.eff 
a   0.07904428 0.03935658 0.02335975 0.05045 0.072255 0.10020 0.1738025 1.001631 2900 
b   0.03557510 0.01966588 0.00946790 0.02168 0.031765 0.04493 0.0842925 1.001946 2000 
c   7.05248063 2.62808107 3.24794998 5.19200 6.603000 8.43600 13.3402499 1.001574 5600 
deviance 61.80275800 2.58984760 58.86000000 59.93000 61.130000 62.98000 68.6200000 1.006417 720 

, , 2 

       mean   sd  2.5%  25%  50%  75%  97.5%  Rhat n.eff 
a   0.1992745 0.08962024 0.06359975 0.1340 0.18575 0.2506 0.413400 1.002039 1800 
b   0.1648684 0.07519227 0.05357925 0.1103 0.15280 0.2070 0.341905 1.001840 2200 
c   7.2121513 1.92005357 3.90300000 5.8610 7.05700 8.4000 11.320250 1.001200 8000 
deviance 47.8683650 2.54986184 44.99000000 46.0300 47.22000 48.9900 54.530000 1.001291 5800 

はパラメータa、b、cのための手段を抽出するが方法ですか?私は、次のようなシミュレーション関数でforループを使用することを考えました:

​​3210

これはうまくいかないようです。どのように私がこの情報を抽出することができるかについてのアイデアは、非常に役立つだろう。ありがとう!

+0

結果の 'dput()'を提供してください。 – Chrisss

+0

@Chrisss私は実際に自分の疑問に対する解決策を見つけました。私は下に自分のソリューションを掲載しました。 –

答えて

0

私は実際に自分の質問を解決しました。反復結果の各値は、out [some number]で抽出できます。たとえば、:

out[1]=0.07904428, out[2]=0.03557510, out[3]=7.05248063 , out[4]=61.80275800, out[5]=0.03935658 

したがって、「out」はすべてのデータを「スネア」するようなものです。私の最初の反復の結果を「蛇行」させたら、それは私の2番目の反復に移動しました。各反復は常に同じ次元であるため、抽出したい各値は常に出発点から36ヵ所離れています。この実現によって、私が望む値だけを得ることができるfor-loopを作成することが可能になりました。

forループ
#This is for acalc 
for (i in 1:length(out)){acalc[i]<-out[(i-1)*36+1]} 
acalc <- as.numeric(na.omit(acalc)) 

#This is for bcalc 
for (i in 1:length(out)){bcalc[i]<-out[(i-1)*36+2]} 
bcalc <- na.omit(bcalc) 

#This is for ccalc 
for (i in 1:length(out)){ccalc[i]<-out[(i-1)*36+3]} 
ccalc <- na.omit(ccalc) 

アウト結果の長さに相当したベクターが得られ、それは唯一各反復におけるパラメータに対応する値が移入されました。これが私がna.omitを使ってNAを取り除いた理由です。

私の質問を審査する時間がかかった人に感謝します。

関連する問題