Search code examples
rjags

finding parameters with jags in r


I'm given the data that contains 1000 numbers from uniform distribution ~[a,b] and I need to use jags in r to find it/ I tried this code

library(arm)
library('rjags')
library(coda)
library(readr)

x <- read.csv(file='C:/Users/Amir/Desktop/אבנר/data analysis/תרגילים/תרגיל 
3/Unif.csv', header=FALSE)
N<-length(x)
y <- x[1:1000,1]

dat<- list("y" = y, "N" = N)
jags.inits <- function() {list (a = -3, b = 30)}
parameters<- c("a", "b")

reg.jags <- jags.model(file='1.txt', data=dat, n.chains = 4, n.adapt = 1000)

update(jags, n.iter=1000)

regression.sim<-coda.samples(reg.jags, variable.names=parameters, 
n.iter=15000)

summary(regression.sim)

and the model is

model {
    for (i in 1:N) {
        y[i] ~ dunif(a, b)
    }

      a ~ dnorm(-5, .0001)
      b ~ dnorm(15, .0001) 

}

but the result are very bad, instead of around [-3,23] I get around [-42,65]

any help?


Solution

  • I can't replicate your findings as your code is not reproducible: we don't have access to your data. But using simulated data with an essentially identical model gives results that seem to be very sensible:

    library('runjags')
    
    m <- 'model{
    
        for(i in 1:N){
            Obs[i] ~ dunif(a, b)
        }
    
        a ~ dnorm(0, 10^-6)
        b ~ dnorm(0, 10^-6)
    
        #data# N, Obs
        #monitor# a, b
        #inits# a, b
    
    }'
    
    N <- 1000
    Obs <- runif(N, 1, 7)
    a <- 0
    b <- 10
    
    results <- run.jags(m)
    plot(results)
    results
    range(Obs)
    

    The upper and lower 95% confidence intervals for a and b respectively are very close to the range of Obs, with mode estimates that are closer to the respective upper and lower 95% CI (the maximum likelihood solution would be exactly at the range of the data), and varying the size of N gives narrower/wider 95% CI. So all as expected.

    My best guess at your problem is that your y are somehow all (or nearly all) missing. If that does not help solve your problem I think you will need to post your dataset.