fixes

pbs-assess · Jan 6, 2024 · 7008acf · 7008acf
1 parent e6a2c25
commit 7008acf
Show file tree

Hide file tree

Showing 4 changed files with 10 additions and 4 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -38,6 +38,7 @@ export(dharma_residuals)
 export(extract_mcmc)
 export(gamma_mix)
 export(gather_sims)
+export(gengamma)
 export(get_cog)
 export(get_crs)
 export(get_index)

diff --git a/R/enum.R b/R/enum.R
@@ -16,7 +16,8 @@
   censored_poisson  = 12,
   gamma_mix = 13,
   lognormal_mix = 14,
-  nbinom2_mix = 15
+  nbinom2_mix = 15,
+  gengamma = 16
 )
 .valid_link <- c(
   identity = 0,

diff --git a/man/families.Rd b/man/families.Rd
diff --git a/src/utils.h b/src/utils.h
@@ -34,9 +34,9 @@ Type dgengamma( Type x,
   // from https://github.com/chjackson/flexsurv-dev/blob/master/src/gengamma.h#L54-L56
   Type y = log(x);
   Type w = (y - mu) / sigma;
-  Type qi = 1/square(Q);
-  Type qw = Q * w;
-  Type logres = -log(sigma*x) + 0.5*log(square(lambda)) * (1 - 2 * qi) + qi * (qw - exp(qw)) - lgamma(qi);
+  Type qi = pow(Q, -2);
+  Type qw = Q * w;                 // 0.5*log(pow(x,2)) as trick for abs(log(x))
+  Type logres = -log(sigma*x) + 0.5*log(pow(lambda,2)) * (1 - 2 * qi) + qi * (qw - exp(qw)) - lgamma(qi);
 
   // return stuff
   if(give_log) return logres; else return exp(logres);