partial work on the Raftery model

Author: maxbiostat

code/Raftery_1988.r

@@ -0,0 +1,18 @@
+library(cmdstanr)
+library(rstan)
+stanfit <- function(fit) rstan::read_stan_csv(fit$output_files())
+
+
+impala <- c(15, 20, 21, 23, 26)
+waterbuck <- c(53, 57, 66, 67, 72)
+
+stan.data <- list(x = impala, n = length(impala), kappa_1 = 1, kappa_2 = 1)
+
+raftery  <- cmdstanr::cmdstan_model("stan/raftery.stan")
+
+fit <- stanfit(raftery$sample(data = stan.data,
+                                  max_treedepth = 15,
+                                  adapt_delta = .99,
+                                  parallel_chains = 4))
+
+fit

code/stan/raftery.stan

@@ -0,0 +1,22 @@
+data{
+  int<lower=1> n;
+  real<lower=0> x[n];
+  real<lower=0> kappa_1;
+  real<lower=0> kappa_2;
+}
+transformed data{
+  real S = sum(x);
+}
+parameters{
+  real<lower=max(x)> N;
+  real<lower=0, upper=1> theta;
+}
+transformed parameters{
+  real lPbar = 0.0;
+  for(j in 1:n) lPbar += lchoose(N, x[j]);
+}
+model{
+ target += -lgamma(N + 1) + lgamma(N + kappa_1) + lPbar;
+ target += (-N + S)*log(theta) + (n*N -S)*log1m(theta)
+ -(N + kappa_1)*log(1/theta + kappa_2);
+}

code/stan/raftery_wrong.stan

@@ -0,0 +1,25 @@
+functions{
+  real bin_lpdf(real x, real n, real p){
+    real ans  = lchoose(n, x) + x*log(p) + (n-x)*log1m(p);
+    return(ans);
+  }
+  real pois_lpdf(real z, real mu){
+    real ans = -mu + z*log(mu) - lgamma(z + 1);
+    return(ans);
+  }
+}
+data{
+  int<lower=1> n;
+  real<lower=0> x[n];
+  real<lower=0> kappa_1;
+  real<lower=0> kappa_2;
+}
+parameters{
+  real<lower=max(x)> N;
+  real<lower=0, upper=1> theta;
+  real<lower=0> mu;
+}
+model{
+  for(j in 1:n) target += bin_lpdf(x[j] | N, theta);
+  target += pois_lpdf(N | mu);
+}