Thanks for helping thanoon.
I used the if() {
....
} else {
.....
}
and it worked better
On 9 February 2016 at 15:56, thanoon younis <thanoon.youni...@gmail.com>
wrote:
> Hi,
>
> Try to check the function only
> if (runif(1) <= accept.prob) beta.cand else beta.curr
> something is missing here. Why you take only runif(1).
>
> Regards
>
> On 9 February 2016 at 16:01, Maram SAlem <marammagdysa...@gmail.com>
> wrote:
>
>> Hi all,
>>
>> I'm trying to write a function to implement a Metropolis-within-Gibbs
>> algorithm for two parameters.I'm including a naive version here so as to
>> be
>> able to spot the error I got. So I first generate the vectors, X and R,
>> that will help to start the algorithm using (for example):
>>
>> n=8; m=5; p=0.1; t=0.9 ; JH=10;
>>
>> R <- numeric(m)
>>
>> W <- numeric(m)
>>
>> V <- numeric(m)
>>
>> U <- numeric(m)
>>
>> X <- numeric(m)
>>
>> Bay.alpha<- numeric (JH)
>>
>> Bay.beta<- numeric (JH)
>>
>> Bay.Surv <- numeric (JH)
>>
>> hyp=c(3,15,6,22.5)
>>
>> theta<-c(0.2,2)
>>
>>
>> alpha.curr<-theta[1]
>>
>> beta.curr<- theta[2]
>>
>>
>>
>> R[1]<-rbinom(1, n-m, p)
>>
>> for (i in 2:m-1) {
>>
>> R[i]<-rbinom(1,n-m-sum(R[1:i-1]),p)
>>
>> }
>>
>> R[m]<-n-m-sum(R[1:m-1])
>>
>> W<-runif(m, min = 0, max = 1)
>>
>> for (i in 1:m){
>>
>> V[i]<-W[i]^(1/(i+sum(R[(m-i+1):m])))
>>
>> }
>>
>> for (i in 1:m){
>>
>> U[i]<- 1- prod(V[(m-i+1):m])
>>
>> }
>>
>> for (i in 1:m){
>>
>> X[i]<- ((-1/theta[1])*log(1-U[i]))^(1/theta[2])
>>
>> }
>>
>>
>> Then, I defined three functions 1- alpha.update() for updating alpha
>> (Gibbs
>> step) 2- bettarg(), for the target distribution of beta 3- beta.update()
>> for updating beta using a Metropolis Hastings technique.
>>
>>
>> alpha.update=function(X, R, alpha.curr, beta.curr ,m, hyp)
>>
>> {
>>
>> o<-numeric(m)
>>
>> for (i in 1:m) {
>>
>> o[i]<- (1+R[i])*((X[i])^(beta.curr))
>>
>> }
>>
>> sh<-sum(o) + hyp[2] + (hyp[4]* beta.curr)
>>
>> rg<-rgamma(1, shape= m+hyp[1]+hyp[3] , rate = sh )
>>
>> return(rg)
>>
>> }
>>
>>
>>
>> bettarg<- function(X, R, alpha.curr, beta.curr ,m, hyp)
>>
>> {
>>
>> o<-numeric(m)
>>
>> for (i in 1:m) {
>>
>> o[i]<- (1+R[i])*((X[i])^( beta.curr))
>>
>> }
>>
>> bt<- beta.curr ^(m+hyp[3]-1) * prod((X)^( beta.curr -1)) *
>> exp(-1*alpha.curr *(sum(o) + (hyp[4]* beta.curr)))
>>
>> return(bt)
>>
>> }
>>
>> beta.update<- function (X, R, alpha.curr, beta.curr ,m, hyp,
>> cand.sd)
>>
>> {
>>
>> beta.cand<- rnorm(1, mean = beta.curr, sd = cand.sd)
>>
>> AB<- bettarg (X, R, alpha.curr, beta.curr = beta.cand ,m, hyp)
>>
>> CD<- bettarg (X, R, alpha.curr, beta.curr ,m, hyp)
>>
>> accept.prob <- AB/CD
>>
>> if (runif(1) <= accept.prob)
>>
>> beta.cand
>>
>> else beta.curr
>>
>> }
>>
>>
>> Then I started a tiny chain using ten iterations but got this ERROR:
>>
>>
>> for (k in 1:JH)
>>
>> {
>>
>> alpha.curr<- alpha.update(X, R, alpha.curr, beta.curr ,m, hyp)
>>
>> Bay.alpha[k]<-alpha.curr
>>
>> beta.curr<- beta.update (X, R, alpha.curr, beta.curr ,m, hyp, cand.sd=2)
>>
>> Bay.beta[k]<- beta.curr
>>
>> Bay.Surv [k]<- exp (-1*Bay.alpha[k] * (t)^Bay.beta[k])
>>
>> }
>>
>>
>> Error in if (runif(1) <= accept.prob) beta.cand else beta.curr :
>>
>> missing value where TRUE/FALSE needed
>>
>> In addition: Warning message:
>>
>> In rgamma(1, shape = m + hyp[1] + hyp[3], rate = sh) : NAs produced
>>
>>
>> I ran the code several times for only one iteration but everything was
>> fine
>> with no errors or warnings, so I don't know from where does the missing
>> value/ NA come from?
>>
>> In addition, I want to calculate the acceptance rate for the
>> Metropolis-Hastings step, is this possible?
>>
>> Any help would be appreciated.
>> Thanks,
>>
>> Maram Salem
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.
>>
>
>
>
> --
> Thanoon Y. Thanoon
> PhD Candidate
> Department of Mathematical Sciences
> Faculty of Science
> University Technology Malaysia, UTM
> E.Mail: thanoon.youni...@gmail.com
> E.Mail: dawn_praye...@yahoo.com
> Facebook:Thanoon Younis AL-Shakerchy
> Twitter: Thanoon Alshakerchy
> H.P:00601127550205
>
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.
