require(rstan)
x=c(3,1,2,5); y=unname(table(x))
dat=list(x=x, y=y, Ndata=length(x))
####
model <- '
data {
int<lower=0> Ndata; int<lower=0> y[Ndata]; int<lower=0> x[Ndata];
}
parameters {
real<lower=0, upper=1> theta;
}
model {
for(i in 1:Ndata){
y[i] ~ poisson(Ndata*(theta^x[i] - theta^(x[i]+1)));
}
}
'
####
fit <- stan(model_code = model, data = dat, iter = 100000, chains = 6, warmup=500)
fit2<-extract(fit)
hist(fit2$theta, breaks=seq(0,1,by=.01))
mean(fit2$theta)
median(fit2$theta)