Skip to content

Commit

Permalink
remove some redundant nmixture code
Browse files Browse the repository at this point in the history
  • Loading branch information
Nicholas Clark committed Oct 31, 2024
1 parent fc765b0 commit 73ab608
Show file tree
Hide file tree
Showing 4 changed files with 0 additions and 87 deletions.
87 changes: 0 additions & 87 deletions R/add_nmixture.R
Original file line number Diff line number Diff line change
Expand Up @@ -114,10 +114,6 @@ add_nmixture = function(model_file,
"Y_max[k] = max(flat_ys[K_inds[k, K_starts[k] : K_stops[k]]]);\n",
"N_max[k] = max(cap[K_inds[k, K_starts[k] : K_stops[k]]]);\n",
"}")
} else {
model_file[grep("transformed data {", model_file, fixed = TRUE)] <-
paste0("transformed data {\n",
"matrix[total_obs, num_basis] X_ordered = X[ytimes_array, : ];")
}
model_file <- readLines(textConnection(model_file), n = -1)

Expand Down Expand Up @@ -379,26 +375,6 @@ add_nmix_model = function(model_file,
model_file, fixed = TRUE)]
model_file <- readLines(textConnection(model_file), n = -1)

} else {
model_file[grep('vector[n_nonmissing] flat_trends;',
model_file, fixed = TRUE)] <- paste0('vector[n_nonmissing] flat_trends;\n',
'vector[n_nonmissing] flat_ps;\n',
'int flat_caps[n_nonmissing];')
model_file <- readLines(textConnection(model_file), n = -1)

model_file[grep('flat_trends = (to_vector(trend))[obs_ind];',
model_file, fixed = TRUE)] <- paste0('flat_trends = (to_vector(trend))[obs_ind];\n',
'flat_ps = p[obs_ind];\n',
'flat_caps = cap[obs_ind];')
model_file <- readLines(textConnection(model_file), n = -1)

model_file[grep('flat_ys ~ poisson_log_glm(append_col(flat_xs, flat_trends),',
model_file, fixed = TRUE)] <- paste0('for (i in 1:n_nonmissing){\n',
'target += pb_lpmf(flat_ys[i] | flat_caps[i], flat_trends[i], flat_ps[i]);\n',
'}')
model_file <- model_file[-grep('0.0,append_row(b, 1.0));',
model_file, fixed = TRUE)]
model_file <- readLines(textConnection(model_file), n = -1)
}

return(model_file)
Expand All @@ -407,69 +383,6 @@ add_nmix_model = function(model_file,
add_nmix_functions = function(model_file,
trend_map,
nmix_trendmap){
if(nmix_trendmap){
# If trend_map supplied, no modifications needed
} else {
if(any(grepl('functions {', model_file, fixed = TRUE))){
model_file[grep('functions {', model_file, fixed = TRUE)] <-
paste0('functions {\n',
'/* Functions to return the log probability of a Poisson Binomial Mixture */\n',
'/* see Bollen et al 2023 for details (https://doi.org/10.1002/ece3.10595)*/\n',
'real poisbin_lpmf(int count, int k, real lambda, real p) {\n',
'if (count > k) {\n',
'return negative_infinity();\n',
'}\n',
'return poisson_log_lpmf(k | lambda) + binomial_logit_lupmf(count | k, p);\n',
'}\n',
'vector pb_logp(int count, int max_k,\n',
'real lambda, real p) {\n',
'if (max_k < count)\n',
'reject("cap variable max_k must be >= observed counts");\n',
'vector[max_k + 1] lp;\n',
'for (k in 0:(count - 1))\n',
'lp[k + 1] = negative_infinity();\n',
'for (k in count:max_k)\n',
'lp[k + 1] = poisbin_lpmf(count | k, lambda, p);\n',
'return lp;\n',
'}\n',
'real pb_lpmf(int count, int max_k,\n',
'real lambda, real p) {\n',
'vector[max_k + 1] lp;\n',
'lp = pb_logp(count, max_k, lambda, p);\n',
'return log_sum_exp(lp);\n',
'}')
} else {
model_file[grep('Stan model code', model_file)] <-
paste0('// Stan model code generated by package mvgam\n',
'functions {\n',
'/* Functions to return the log probability of a Poisson Binomial Mixture */\n',
'/* see Bollen et al 2023 for details (https://doi.org/10.1002/ece3.10595)*/\n',
'real poisbin_lpmf(int count, int k, real lambda, real p) {\n',
'if (count > k) {\n',
'return negative_infinity();\n',
'}\n',
'return poisson_log_lpmf(k | lambda) + binomial_logit_lupmf(count | k, p);\n',
'}\n',
'vector pb_logp(int count, int max_k,\n',
'real lambda, real p) {\n',
'if (max_k < count)\n',
'reject("cap variable max_k must be >= observed counts");\n',
'vector[max_k + 1] lp;\n',
'for (k in 0:(count - 1))\n',
'lp[k + 1] = negative_infinity();\n',
'for (k in count:max_k)\n',
'lp[k + 1] = poisbin_lpmf(count | k, lambda, p);\n',
'return lp;\n',
'}\n',
'real pb_lpmf(int count, int max_k,\n',
'real lambda, real p) {\n',
'vector[max_k + 1] lp;\n',
'lp = pb_logp(count, max_k, lambda, p);\n',
'return log_sum_exp(lp);\n',
'}\n}\n')
}
}

model_file <- readLines(textConnection(model_file), n = -1)
return(model_file)
}
Expand Down
Binary file added src/RcppExports.o
Binary file not shown.
Binary file added src/trend_funs.o
Binary file not shown.
Binary file modified tests/testthat/Rplots.pdf
Binary file not shown.

0 comments on commit 73ab608

Please sign in to comment.