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