initial values for a proportional estimate

MDP
2013-05-09
2013-05-10
  • MDP
    MDP
    2013-05-09

    I am working through a basic problem from Kruschke's book in preparation to modify the code for my particular problem Specifically, I am trying to reproduce the binomial proportion example as follows:

    Given 14 coin flips, 11 flips are "heads" or successful. Estimate the proportion of successes from the data below:

    coin_data <- c(r(1,11), r(0,3))
    N <- length(coin_data)
    binom_estimate <- "model{
    
        for(i in 1:N){
            y[i] ~ dbern(theta)
        }
    # Priors
        theta ~ dbin(priorA, priorB)
        priorA <- 1
        priorB <- 1
    }
    "
    writeLines(binom_estimate, con="binom_estimate.txt")
    

    And, here is the implementation... cutting out some of the stuff Kruschke does for
    simplification.

    # required packages
    library(rjags)
    
    # set up the bayesian estimate
    data_list <- list(y = coin_data,
                      N = N)
    
    # run the model
    bayes_estimate <- jags.model("binom_estimate.txt",
                                 data = data_list,
                                 n.chains = 3, 
                                 n.adapt = 500)
    
    update(bayes_estimate)
    
    bayes_estimated_results <- coda.samples(bayes_estimate,
                                            c("theta", "priorA", "priorB"),
                                            n.iter=1000)
    
    summary(bayes_estimated_results)
    

    My result, with my current version of R2jags = 0.03-08 and JAGS = 3.3.0 is:

    Error in node y[12]
    Observed node inconsistent with unobserved parents at initialization
    

    If anyone has any suggestion on how I can fix this or how to provide useful initializations,
    I would love to read them. I have read through the threads posted on this type of topic and
    tried the solutions that seem to fit my particular problem, but haven't been able to get it
    work.

    Best,

    MDP

     
  • Martyn Plummer
    Martyn Plummer
    2013-05-09

    There are a couple of typos in your code. Firt, there is no r function in R, but there is a rep function, which I'm guessing is what you meant

    Secondly, the line

    theta ~ dbin(priorA, priorB)
    

    in your BUGS code should be replaced by

    theta ~ dbeta(priorA, priorB)
    

    Then everything should work fine.

     
  • MDP
    MDP
    2013-05-10

    UGH!!!! I feel like such a heal...THANK YOU.