1

I am attempting to compare two beta distributions under a Bayesian framework. (My data are survival rates, so they fall between 0 and 1 and are best fit with a beta distribution. After estimating alpha and beta, I will moment match to compare means.) However, the beta likelihood is returning "0" for the given initial values, which makes estimating the accept/reject ratio for MCMC in the given code (at the end of my question below) impossible to calculate.

I am gathering a rudimentary understanding of why this occurs, but am in a bit of a hurry and would rather not dig too far into the technical literature, so am hoping someone might have an easily-implementable solution I overlooked. One suggestion would be to somehow incorporate code similar to that which I found in a "Bayesian Model for Practicing Ecologists" lab (which did not include solutions):

L.vector = dbeta(y,a,b, log=TRUE) like=sum(L.vector[is.finite(L.vector)])

However, I am not sure where to insert this code, whether using "log=TRUE" is appropriate here, or whether there would even be enough instances with finite values to succeed.

My code is below, thank you SO MUCH to anyone who takes the time to give suggestions.

likelihood <- function(parameters){
  alpha.1=parameters[1]; beta.1=parameters[2]; alpha.2=parameters[3]; beta.2=parameters[4]
  prod(dbeta(data.1, alpha.1, beta.1)) * prod(dbeta(data.2, alpha.2, beta.2))
}

prior <- function(parameters){
  alpha.1=parameters[1]; beta.1=parameters[2]; alpha.2 =parameters[3]; beta.2=parameters[4]
  dgamma(alpha.1, shape.prior, rate.prior) * dgamma(beta.1, shape.prior, rate.prior) * dgamma(alpha.2, shape.prior, rate.prior) * dgamma(beta.2, shape.prior, rate.prior)}

posterior <- function(parameters) {likelihood(parameters) * prior(parameters)}

#starting values
alpha.1 = .5; beta.1 = .5; alpha.2 = .5; beta.2 = .5
parameters <- c(alpha.1, beta.1, alpha.2, beta.2)

#MCMC w/Metropolis
n.iter <- 100000
results <- matrix(0, nrow=n.iter, ncol=4)
results[1, ] <- parameters
for (iteration in 2:n.iter){
  candidate <- parameters + runif(4)
  ratio <- posterior(candidate)/posterior(parameters)
  if (runif(1,0,1) < ratio) {
    parameters <- candidate 
  results[iteration, ] <- parameters
  }
}

Also, thanks to Kruschke 2013 and the person who gave code at Bayesian equivalent of two sample t-test? for getting me to this point!

user6776473
  • 19
  • 3
  • Hi Xi'an, thanks for that suggestion! I suppose I could try logging everything, then I assume it would work to back-transform the resulting parameters to get alpha and beta? Am I thinking correctly? – user6776473 Aug 10 '20 at 03:32

1 Answers1

2

To check whether using a logarithmic computation of all functions (not of parameters) the lines you need change in your R code are those

sum(dbeta(data.1, alpha.1, beta.1,log=TRUE)) + 
sum(dbeta(data.2, alpha.2, beta.2,log=TRUE))

then

    dgamma(alpha.1, shape.prior, rate.prior,log=TRUE) + 
    dgamma(beta.1, shape.prior, rate.prior,log=TRUE)  + 
    dgamma(alpha.2, shape.prior, rate.prior,log=TRUE) +
    dgamma(beta.2, shape.prior, rate.prior,log=TRUE)}

then

    posterior <- function(parameters) {likelihood(parameters) + prior(parameters)}

and

    ratio <- posterior(candidate)-posterior(parameters)
    if (log(runif(1,0,1)) < ratio) {

Furthermore, your Metropolis-Hastings step is not correct, it should be

    if (log(runif(1,0,1)) < ratio) parameters <- candidate
     results[iteration, ] <- parameters

so that the rejection of the candidate induces a repetition of the current value of the parameters.

And the proposal is also incorrect

    candidate <- parameters + runif(4)

should be

    candidate <- parameters + runif(4,-1,1)

for the Markov chain to be reversible.

Xi'an
  • 90,397
  • 9
  • 157
  • 575
  • Thank you Xi'an! I tried all of this, but am still receiving this error: Error in if (log(runif(1, 0, 1)) < ratio) parameters – user6776473 Aug 11 '20 at 18:58
  • You need to reject every negative proposal before computing the ratio: ` if (min(candidate)>0){` – Xi'an Aug 12 '20 at 06:20
  • See [full code](https://shorturl.at/hjAB5) that runs with no error. – Xi'an Aug 12 '20 at 06:24
  • Xi'an, I cannot thank you enough! This works. Many of the parameters in the results matrix come out at 0 though; there are many rows of 0's, which seems odd. In addition, it could be a virus on my computer, but the link to your "full code" takes me to a short URL page so I cannot view it. I inserted that new line right before the ratio bit as suggested, and it ran, so I am assuming that is correct. – user6776473 Aug 16 '20 at 00:06