Hi all, I appreciate the help this list has given me before. I have a
question which has been perplexing me. I have been working on doing a
Bayesian calculating inserting studies sequentially after using a
non-informative prior to get a meta-analysis type result. I created a
function using three iterations of this, my code is below. I insert prior
mean and precision (I add precision manually because I want it to be zero
but if I include zero as prior variance I get an error cant divide by zero.
Now my question is instead of having the three iterations below, is there
anyway to create a loop, the only problem is I want to have elements from
the previous posterior to be the new prior and now I cant figure out how to
do the code below in a loop. The data below is dummy data, I used a starting
mu of 1, and starting precision of 0.
bayes.analysis.treat<-function(mu0,p0){
n1 = 5
n2 = 10
n3 = 15
ybar1 = 12
ybar2 = 13
ybar3 = 14
sd1 = 2
sd2 = 3
sd3 = 4
#posterior 1
var1 = sd1^2 #sample variance
p1 = n1/var1 #sample precision
p1n = p0+p1
mu1 = ((p0)/(p1n)*mu0)+((p1)/(p1n))*ybar1
sigma1 = 1/sqrt(p1n)
#posterior 2
var2 = sd2^2 #sample variance
p2 = n2/var2 #sample precision
p2n = p1n+p2
mu2 = ((p1n)/(p2n)*mu1)+((p2)/(p2n))*ybar2
sigma2 = 1/sqrt(p2n)
#posterior 3
var3 = sd3^2 #sample variance
p3 = n3/var3 #sample precision
p3n = p2n+p3
mu3 = ((p2n)/(p3n)*mu2)+((p3)/(p3n))*ybar3
sigma3 = 1/sqrt(p3n)
posterior<-cbind(rbind(mu1,mu2,mu3),rbind(sigma1,sigma2,sigma3))
rownames(posterior)<-c("Posterior 1", "Posterior 2",
"Posterior 3")
colnames(posterior)<-c("Mu", "SD")
return(posterior)}
-------------------------------------------
Joe King, M.A.
Ph.D. Student
University of Washington - Seattle
206-913-2912
<mailto:jp@joepking.com> jp@joepking.com
-------------------------------------------
"Never throughout history has a man who lived a life of ease left a name
worth remembering." --Theodore Roosevelt
[[alternative HTML version deleted]]
On Jul 17, 2010, at 9:09 PM, Joe P King wrote:> Hi all, I appreciate the help this list has given me before. I have a > question which has been perplexing me. I have been working on doing a > Bayesian calculating inserting studies sequentially after using a > non-informative prior to get a meta-analysis type result. I created a > function using three iterations of this, my code is below. I insert > prior > mean and precision (I add precision manually because I want it to be > zero > but if I include zero as prior variance I get an error cant divide > by zero. > Now my question is instead of having the three iterations below, is > there > anyway to create a loop, the only problem is I want to have elements > from > the previous posterior to be the new prior and now I cant figure out > how to > do the code below in a loop.Perhaps if you set the problem up with vectors, it would become more obvious how to do it with indexing or loops: n <- c(5, 10, 15) ybar <- c(12,12,14) # and so on... -- David> The data below is dummy data, I used a starting > mu of 1, and starting precision of 0. > > > > bayes.analysis.treat<-function(mu0,p0){ > > n1 = 5 > > n2 = 10 > > n3 = 15 > > ybar1 = 12 > > ybar2 = 13 > > ybar3 = 14 > > sd1 = 2 > > sd2 = 3 > > sd3 = 4 > > #posterior 1 > > var1 = sd1^2 #sample variance > > p1 = n1/var1 #sample precision > > p1n = p0+p1 > > mu1 = ((p0)/(p1n)*mu0)+((p1)/(p1n))*ybar1 > > sigma1 = 1/sqrt(p1n) > > #posterior 2 > > var2 = sd2^2 #sample variance > > p2 = n2/var2 #sample precision > > p2n = p1n+p2 > > mu2 = ((p1n)/(p2n)*mu1)+((p2)/(p2n))*ybar2 > > sigma2 = 1/sqrt(p2n) > > #posterior 3 > > var3 = sd3^2 #sample variance > > p3 = n3/var3 #sample precision > > p3n = p2n+p3 > > mu3 = ((p2n)/(p3n)*mu2)+((p3)/(p3n))*ybar3 > > sigma3 = 1/sqrt(p3n) > > posterior<-cbind(rbind(mu1,mu2,mu3),rbind(sigma1,sigma2,sigma3)) > > rownames(posterior)<-c("Posterior 1", "Posterior 2", "Posterior 3") > > colnames(posterior)<-c("Mu", "SD") > > return(posterior)} > > > > ------------------------------------------- > > Joe King, M.A. > > Ph.D. Student > > University of Washington - Seattle > > 206-913-2912 > > <mailto:jp at joepking.com> jp at joepking.com > > ------------------------------------------- > > "Never throughout history has a man who lived a life of ease left a > name > worth remembering." --Theodore Roosevelt > > > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.David Winsemius, MD West Hartford, CT
I tried that, this is what I tried, but I only get it to do one iteration
and then it wont cycle back. I am not sure how to tell it how to look at the
previous precision to add to the current new sample precision to get the new
precision. This is my first try at a loop so I am not sure if I am doing
anything right, sorry about the elementary question.
n=c(5,10,15)
ybar=c(12,13,14)
sd=c(2,3,4)
mutotals<-matrix(nrow=length(n), ncol=1)
mu=1;p0=0
for (i in 1:length(n)){
var = sd[i]^2 #sample variance
p = n[i]/var #sample precision
p0[i]= i
pn = p0[i]+p[i]
mu.out = ((p0[i])/(pn)*mu[i])+((p)/(pn))*ybar[i]
sigma[i] = 1/sqrt(pn[i])
mutotals[i,]<-mu.out[i]
}
Joe King
206-913-2912
jp at joepking.com
"Never throughout history has a man who lived a life of ease left a name
worth remembering." --Theodore Roosevelt
-----Original Message-----
From: David Winsemius [mailto:dwinsemius at comcast.net]
Sent: Sunday, July 18, 2010 5:36 AM
To: Joe P King
Cc: r-help at r-project.org
Subject: Re: [R] loop troubles
On Jul 17, 2010, at 9:09 PM, Joe P King wrote:
> Hi all, I appreciate the help this list has given me before. I have a
> question which has been perplexing me. I have been working on doing a
> Bayesian calculating inserting studies sequentially after using a
> non-informative prior to get a meta-analysis type result. I created a
> function using three iterations of this, my code is below. I insert
> prior
> mean and precision (I add precision manually because I want it to be
> zero
> but if I include zero as prior variance I get an error cant divide
> by zero.
> Now my question is instead of having the three iterations below, is
> there
> anyway to create a loop, the only problem is I want to have elements
> from
> the previous posterior to be the new prior and now I cant figure out
> how to
> do the code below in a loop.
Perhaps if you set the problem up with vectors, it would become more
obvious how to do it with indexing or loops:
n <- c(5, 10, 15)
ybar <- c(12,12,14) # and so on...
--
David> The data below is dummy data, I used a starting
> mu of 1, and starting precision of 0.
>
>
>
> bayes.analysis.treat<-function(mu0,p0){
>
> n1 = 5
>
> n2 = 10
>
> n3 = 15
>
> ybar1 = 12
>
> ybar2 = 13
>
> ybar3 = 14
>
> sd1 = 2
>
> sd2 = 3
>
> sd3 = 4
>
> #posterior 1
>
> var1 = sd1^2 #sample variance
>
> p1 = n1/var1 #sample precision
>
> p1n = p0+p1
>
> mu1 = ((p0)/(p1n)*mu0)+((p1)/(p1n))*ybar1
>
> sigma1 = 1/sqrt(p1n)
>
> #posterior 2
>
> var2 = sd2^2 #sample variance
>
> p2 = n2/var2 #sample precision
>
> p2n = p1n+p2
>
> mu2 = ((p1n)/(p2n)*mu1)+((p2)/(p2n))*ybar2
>
> sigma2 = 1/sqrt(p2n)
>
> #posterior 3
>
> var3 = sd3^2 #sample variance
>
> p3 = n3/var3 #sample precision
>
> p3n = p2n+p3
>
> mu3 = ((p2n)/(p3n)*mu2)+((p3)/(p3n))*ybar3
>
> sigma3 = 1/sqrt(p3n)
>
> posterior<-cbind(rbind(mu1,mu2,mu3),rbind(sigma1,sigma2,sigma3))
>
> rownames(posterior)<-c("Posterior 1", "Posterior 2",
"Posterior 3")
>
> colnames(posterior)<-c("Mu", "SD")
>
> return(posterior)}
>
>
>
> -------------------------------------------
>
> Joe King, M.A.
>
> Ph.D. Student
>
> University of Washington - Seattle
>
> 206-913-2912
>
> <mailto:jp at joepking.com> jp at joepking.com
>
> -------------------------------------------
>
> "Never throughout history has a man who lived a life of ease left a
> name
> worth remembering." --Theodore Roosevelt
>
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
West Hartford, CT