Version fixes and minor typos (in progress)

.gitignore

@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

COMP_Stan.Rproj

@@ -0,0 +1,13 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: ISO8859-1
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX

fit_simple_COMP.r

@@ -12,6 +12,7 @@ MaxIter <- 1E4
 nobs <- 1000
 # set.seed(666)
 Y <- COMPoissonReg::rcmp(n = nobs, lambda = Mu, nu = Nu)
+x11(width = 8, height = 6)
 hist(Y, probability = TRUE)
 
 cY <- compress_counts(Y)

fit_simple_COMP_parametrisation_comparison.r

@@ -12,6 +12,7 @@ MaxIter <- 1E4
 nobs <- 1000
 # set.seed(666)
 Y <- COMPoissonReg::rcmp(n = nobs, lambda = Mu, nu = Nu)
+x11(width = 8, height = 6)
 hist(Y, probability = TRUE)
 
 cY <- compress_counts(Y)

stan/ErrorBounding.stan

@@ -1,6 +1,6 @@
 /* Modified brms implementation that sums all terms at once instead of doing it in batches */
 // Taken from https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_com_poisson.stan
-real[] log_Z_COMP_EB(real log_mu, real nu, real eps, int M) {
+array[] real log_Z_COMP_EB(real log_mu, real nu, real eps, int M) {
   int k = 2;
   real leps = log(eps);
   vector[M] log_Z_terms;
@@ -15,6 +15,6 @@ real[] log_Z_COMP_EB(real log_mu, real nu, real eps, int M) {
     reject("nu must be finite");
   }
   return(
-    infiniteErrorBoundingPairs({log_mu, nu}, eps, M, log(0), 0)
+    infiniteErrorBoundingPairs({log_mu, nu}, eps, M, 0)
   );
 }

stan/Fixed.stan

@@ -1,6 +1,6 @@
 /* Modified brms implementation that sums all terms at once instead of doing it in batches */
 // Taken from https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_com_poisson.stan
-real[] log_Z_COMP_fixed(real log_mu, real nu, int Nmax, int M) {
+array[] real log_Z_COMP_fixed(real log_mu, real nu, int Nmax, int M) {
   vector[M] log_Z_terms;
   if (nu == 1) {
     return {exp(log_mu), 0};

stan/Fixed_2.stan

@@ -1,6 +1,6 @@
 /* Modified brms implementation that sums all terms at once instead of doing it in batches */
 // Taken from https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_com_poisson.stan
-real[] log_Z_COMP_2_fixed(real log_mu, real nu, int Nmax, int M) {
+array[] real log_Z_COMP_2_fixed(real log_mu, real nu, int Nmax, int M) {
   vector[M] log_Z_terms;
   if (nu == 1) {
     return {exp(log_mu), 0};

stan/Gaunt.stan

@@ -4,7 +4,7 @@ real log_Z_COMP_Asymp(real log_mu, real nu) {
   // Based on equations (4) and (31) of doi:10.1007/s10463-017-0629-6
   real nu2 = nu^2;
   real log_common = log(nu) + log_mu/nu;
-  real resids[4];
+  array[4] real resids;
   real ans;
   real lcte = (nu * exp(log_mu/nu)) - ( (nu-1)/(2*nu)* log_mu + (nu-1)/2*log(2*pi()) + 0.5 *log(nu));
   real c_1 = (nu2-1)/24;

stan/MCMC_COMP_2_Batches.stan

@@ -4,8 +4,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -21,7 +21,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = infiniteBatches({log_mu, nu}, batchSize, eps, M, 0);
+  array[2] real log_norm_const = infiniteBatches({log_mu, nu}, batchSize, eps, M, 0);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_2_ErrorBounding.stan

@@ -4,8 +4,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -14,16 +14,13 @@ data{
   real<lower=0> eps;
   int<lower=0> M;
 }
-transformed data{
-  real logL = log(0);
-}
 parameters{
   real<lower=0> mu;
   real<lower=0> nu;
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = infiniteErrorBoundingPairs({log_mu, nu}, eps, M, logL, 0);
+  array[2] real log_norm_const = infiniteErrorBoundingPairs({log_mu, nu}, eps, M, 0);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_2_Threshold.stan

@@ -5,8 +5,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -21,7 +21,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = log_Z_COMP_2_brms(log_mu, nu, eps, M);
+  array[2] real log_norm_const = log_Z_COMP_2_brms(log_mu, nu, eps, M);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_2_brms.stan

@@ -4,8 +4,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -19,7 +19,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = log_Z_COMP_brms2(log_mu, nu, eps, M);
+  array[2] real log_norm_const = log_Z_COMP_brms2(log_mu, nu, eps, M);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_2_brms_bulk.stan

@@ -4,8 +4,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -20,7 +20,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = log_Z_COMP_brms_bulk2(log_mu, nu, eps, M);
+  array[2] real log_norm_const = log_Z_COMP_brms_bulk2(log_mu, nu, eps, M);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_2_fixedCap.stan

@@ -5,8 +5,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -21,7 +21,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = log_Z_COMP_2_fixed(log_mu, nu, N_cap, M);
+  array[2] real log_norm_const = log_Z_COMP_2_fixed(log_mu, nu, N_cap, M);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_Batches.stan

@@ -4,8 +4,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -20,7 +20,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = infiniteBatches({log_mu, nu}, batchSize, eps, M, 0);
+  array[2] real log_norm_const = infiniteBatches({log_mu, nu}, batchSize, eps, M, 0);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_ErrorBounding.stan

@@ -5,8 +5,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -20,7 +20,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = log_Z_COMP_EB(log_mu, nu, eps, M);
+  array[2] real log_norm_const = log_Z_COMP_EB(log_mu, nu, eps, M);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_Threshold.stan

@@ -5,8 +5,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -20,7 +20,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = log_Z_COMP_Threshold(log_mu, nu, eps, M);
+  array[2] real log_norm_const = log_Z_COMP_Threshold(log_mu, nu, eps, M);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_brms.stan

@@ -4,8 +4,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -19,7 +19,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = log_Z_COMP_brms(log_mu, nu, eps, M);
+  array[2] real log_norm_const = log_Z_COMP_brms(log_mu, nu, eps, M);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_brms_bulk.stan

@@ -4,8 +4,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -19,7 +19,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = log_Z_COMP_brms_bulk(log_mu, nu, eps, M);
+  array[K] real log_norm_const = log_Z_COMP_brms_bulk(log_mu, nu, eps, M);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/MCMC_COMP_fixedCap.stan

@@ -5,8 +5,8 @@ functions{
 }
 data{
   int<lower=0> K;
-  int<lower=0> n[K];
-  int<lower=0> y[K];
+  array[K] int<lower=0> n;
+  array[K] int<lower=0> y;
   int<lower=0> N;
   real<lower=0> s_mu;
   real<lower=0> r_mu;
@@ -20,7 +20,7 @@ parameters{
 }
 transformed parameters{
   real log_mu = log(mu);
-  real log_norm_const[2] = log_Z_COMP_fixed(log_mu, nu, N_cap, M);
+  array[2] real log_norm_const = log_Z_COMP_fixed(log_mu, nu, N_cap, M);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

stan/SumToThreshold.stan

@@ -1,6 +1,6 @@
 /* Modified brms implementation that sums all terms at once instead of doing it in batches */
 // Taken from https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_com_poisson.stan
-real[] log_Z_COMP_Threshold(real log_mu, real nu, real eps, int M) {
+array[] real log_Z_COMP_Threshold(real log_mu, real nu, real eps, int M) {
   int k = 2;
   real leps = log(eps);
   vector[M] log_Z_terms;

stan/SumToThreshold_2.stan

@@ -1,6 +1,6 @@
 /* Modified brms implementation that sums all terms at once instead of doing it in batches */
 // Taken from https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_com_poisson.stan
-real[] log_Z_COMP_2_brms(real log_mu, real nu, real eps, int M) {
+array[] real log_Z_COMP_2_brms(real log_mu, real nu, real eps, int M) {
   int k = 2;
   real leps = log(eps);
   vector[M] log_Z_terms;

stan/brms.stan

@@ -1,6 +1,6 @@
 /* Modified brms implementation with bug fixes and a bit of streamlining. Keeps the 'batch' approach */
 // Taken from https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_com_poisson.stan
-real[] log_Z_COMP_brms(real log_mu, real nu, real eps, int M) {
+array[] real log_Z_COMP_brms(real log_mu, real nu, real eps, int M) {
   real log_Z;
   int k = 2;
   real leps = log(eps);

stan/brms_2.stan

@@ -1,6 +1,6 @@
 /* Modified brms implementation with bug fixes and a bit of streamlining.*/
 // Taken from https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_com_poisson.stan
-real[] log_Z_COMP_brms2(real log_mu, real nu, real eps, int M) {
+array[] real log_Z_COMP_brms2(real log_mu, real nu, real eps, int M) {
   real log_Z;
   int k = 2;
   real leps = log(eps);

stan/brms_bulk.stan

@@ -1,6 +1,6 @@
 /* Modified brms implementation that sums all terms at once instead of doing it in batches */
 // Taken from https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_com_poisson.stan
-real[] log_Z_COMP_brms_bulk(real log_mu, real nu, real eps, int M) {
+array[] real log_Z_COMP_brms_bulk(real log_mu, real nu, real eps, int M) {
   int k = 2;
   real leps = log(eps);
   vector[M] log_Z_terms;

stan/brms_bulk_2.stan

@@ -1,6 +1,6 @@
 /* Modified brms implementation that sums all terms at once instead of doing it in batches */
 // Taken from https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_com_poisson.stan
-real[] log_Z_COMP_brms_bulk2(real log_mu, real nu, real eps, int M) {
+array[] real log_Z_COMP_brms_bulk2(real log_mu, real nu, real eps, int M) {
   int k = 2;
   real leps = log(eps);
   vector[M] log_Z_terms;

stan/comp_2_pmf.stan

@@ -5,6 +5,6 @@ real log_COM_Poisson_2(int k, real log_mu, real nu){
 real COM_Poisson_2_lpmf(int k, real log_mu, real nu, real logZ){
     return nu*(k * log_mu - lgamma(k + 1)) - logZ;
 }
-real logFunction(int n, real[] p){
+real logFunction(int n, array[] real p){
   return(log_COM_Poisson_2(n, p[1], p[2]));
 }

stan/comp_pmf.stan

@@ -5,6 +5,6 @@ real log_COM_Poisson(int k, real log_mu, real nu){
 real COM_Poisson_lpmf(int k, real log_mu, real nu, real logZ){
     return k * log_mu - nu * lgamma(k + 1) - logZ;
 }
-real logFunction(int n, real[] p){
-  return(log_COM_Poisson(n, p[1], p[2]));
+real logFunction(int n, array[] real p){
+  return log_COM_Poisson(n, p[1], p[2]);
 }