2016-10-05 4 views
0

Я относительно новичок в R и выполняю симуляцию MCMC с использованием R2OpenBUGS. Моя модель работает, и я стараюсь делать несколько итераций эффективно. Это мой код для нескольких итераций.Извлечение значений из реплицируемых функций

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 

Есть ли способ, чтобы извлечь средства для параметров а, Ь, с? Я думал об использовании цикла for в функции моделирования, например:

acalc<-numeric(length=iter) 
bcalc<-numeric(length=iter) 
ccalc<-numeric(length=iter) 
for (i in 1:iter){acalc[i]<-out["a","mean"];bcalc[i]<-out["b","mean"];ccalc[i]<-out["c","mean"]} 

Но это не работает. Любые идеи о том, как я мог бы извлечь эту информацию, были бы чрезвычайно полезными. Благодаря!

+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 

Поэтому «выход» - это как один большой вектор, «пронзающий» по всем моим данным. Как только он «просочился» через мои первые результаты итерации, он перешел на вторую итерацию. Поскольку каждая итерация всегда будет иметь одинаковую размерность, каждое значение, которое я хочу извлечь, всегда будет находиться в 36 местах от моей начальной точки. Эта реализация позволила мне создать цикл 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, чтобы избавиться от НС.

Спасибо всем, кто нашел время, чтобы рассмотреть мой вопрос.

Смежные вопросы