CODEWITHSUNDEEP
rjags
נעמי נבו
נעמי נבו

Reputation: 3

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?

Upvotes: 0

Views: 507

Answers (1)

Matt Denwood
Matt Denwood

Reputation: 2583

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.

Upvotes: 1

Related Questions