Extended family=binomialmix

Revised family=binomialmix

gmrflib/GMRFLib.h

@@ -163,7 +163,7 @@ void mkl_dcsrmv(const char *transa, const int *m, const int *k, const double *al
 double cblas_ddoti(const int nz, const double *x, const int *indx, const double *y);
 #endif
 
-#if defined(INLA_WITH_FRAMEWORK_ACCELERATE) 
+#if defined(INLA_WITH_FRAMEWORK_ACCELERATE)
 void vvsqrt(double *, const double *, const int *);
 void vvexp(double *, const double *, const int *);
 void vvexpm1(double *, const double *, const int *);

gmrflib/dot.c

@@ -136,8 +136,17 @@ double GMRFLib_dot_product_serial_mkl_alt(GMRFLib_idxval_tp *__restrict ELM_, do
 
 double GMRFLib_ddot(int n, double *x, double *y)
 {
-	int one = 1;
-	return ddot_(&n, x, &one, y, &one);
+	if (n <= 4L) {
+		double r = 0.0;
+#pragma omp simd reduction(+: r)
+		for (int i = 0; i < n; i++) {
+			r += x[i] * y[i];
+		}
+		return r;
+	} else {
+		int one = 1;
+		return ddot_(&n, x, &one, y, &one);
+	}
 }
 
 

gmrflib/idxval.c

@@ -1318,7 +1318,7 @@ int GMRFLib_str_is_member(GMRFLib_str_tp *hold, char *s, int case_sensitive, int
 		return 0;
 	}
 
-	int (*cmp)(const char *, const char *) =(case_sensitive ? strcmp : strcasecmp);
+	int (*cmp)(const char *, const char *) = (case_sensitive ? strcmp : strcasecmp);
 	for (int i = 0; i < hold->n; i++) {
 		if (cmp(s, hold->str[i]) == 0) {
 			if (idx_match) {

gmrflib/simd.c

@@ -93,7 +93,7 @@ void GMRFLib_exp(int n, double *x, double *y)
 		vdExp(n, x, y);
 #elif defined(INLA_WITH_FRAMEWORK_ACCELERATE)
 		vvexp(y, x, &n);
-#else		
+#else
 		_Pragma("omp simd")
 		    for (int i = 0; i < n; i++) {
 			y[i] = exp(x[i]);

gmrflib/utils.c

@@ -2111,9 +2111,9 @@ int GMRFLib_is_sorted_ddec_plain(int n, double *a)
 
 int GMRFLib_is_sorted(void *a, size_t n, size_t size, int (*cmp)(const void *, const void *))
 {
-	if((cmp ==(void *) GMRFLib_icmp) && size == sizeof(int)) {
+	if ( (cmp == (void *) GMRFLib_icmp) && size == sizeof(int)) {
 		// increasing ints
-		return GMRFLib_is_sorted_iinc(n,(int *) a);
+		return GMRFLib_is_sorted_iinc(n, (int *) a);
 	} else if (cmp == (void *) GMRFLib_icmp_r && size == sizeof(int)) {
 		// decreasing ints
 		return GMRFLib_is_sorted_idec(n, (int *) a);
@@ -2133,15 +2133,15 @@ int GMRFLib_is_sorted(void *a, size_t n, size_t size, int (*cmp)(const void *, c
 void GMRFLib_qsort(void *a, size_t n, size_t size, int (*cmp)(const void *, const void *))
 {
 	// sort if not sorted
-	if(n > 0 && !GMRFLib_is_sorted(a, n, size, cmp)) {
+	if (n > 0 && !GMRFLib_is_sorted(a, n, size, cmp)) {
 		QSORT_FUN(a, n, size, cmp);
 	}
 }
 
 void GMRFLib_qsort2(void *x, size_t nmemb, size_t size_x, void *y, size_t size_y, int (*compar)(const void *, const void *))
 {
-	if(!y) {
-		return (GMRFLib_qsort(x, nmemb, size_x, compar));
+	if (!y) {
+		return(GMRFLib_qsort(x, nmemb, size_x, compar));
 	}
 
 	if (nmemb == 0) {

inlaprog/src/inla-likelihood.c

@@ -597,7 +597,7 @@ int inla_read_data_likelihood(inla_tp *mb, dictionary *UNUSED(ini), int UNUSED(s
 
 	case L_BINOMIALMIX:
 	{
-		assert(ncol_data_all == 13);
+		assert(ncol_data_all == 16);
 		idiv = ncol_data_all;
 		na = ncol_data_all - 2;
 		for (i = 0; i < na; i++) {
@@ -3318,7 +3318,7 @@ int loglikelihood_binomialmix(int thread_id, double *__restrict logll, double *_
 	Data_section_tp *ds = (Data_section_tp *) arg;
 	double y = ds->data_observations.y[idx];
 	double *dat = ds->data_observations.binmix_dat[idx];
-	double n = dat[BINOMIALMIX_NBETA + 2];
+	double n = dat[BINOMIALMIX_NBETA + 4];		       /* z_1, ..., z_11, w_1, w_2 */
 	double ny = n - y;
 
 	if (ISZERO(y) && ISZERO(n)) {
@@ -3338,49 +3338,48 @@ int loglikelihood_binomialmix(int thread_id, double *__restrict logll, double *_
 	double normc = G_norm_const[idx];
 	LINK_INIT;
 
-	int nbeta2 = BINOMIALMIX_NBETA / 2L;		       /* yes, _NBETA is even */
-	double p1_intern = 0.0;
-	for (int i = 0; i < nbeta2; i++) {
-		p1_intern += ds->data_observations.binmix_beta[i][thread_id][0] * dat[i];
-	}
-	double p1 = PREDICTOR_INVERSE_LINK(p1_intern);
-
-	int offset = nbeta2;
-	double p2_intern = 0.0;
-	for (int i = 0; i < nbeta2; i++) {
-		p2_intern += ds->data_observations.binmix_beta[offset + i][thread_id][0] * dat[offset + i];
-	}
-	double p2 = PREDICTOR_INVERSE_LINK(p2_intern);
-
-	double p12 = dat[BINOMIALMIX_NBETA] * p1 + dat[BINOMIALMIX_NBETA + 1] * p2;
-	double w3 = 1.0 - (dat[BINOMIALMIX_NBETA] + dat[BINOMIALMIX_NBETA + 1]);
+	double beta[BINOMIALMIX_NBETA];
+	for (int i = 0; i < BINOMIALMIX_NBETA; i++) {
+		beta[i] = ds->data_observations.binmix_beta[i][thread_id][0];
+	}
+	double beta9 = beta[BINOMIALMIX_NBETA-1];
+
+	int nbeta2 = BINOMIALMIX_NBETA / 2L;
+	double p1_intern = beta9 * dat[BINOMIALMIX_NBETA-1] + GMRFLib_ddot(nbeta2, beta, dat);
+	double p2_intern = beta9 * dat[BINOMIALMIX_NBETA+0] + GMRFLib_ddot(nbeta2, beta + nbeta2, dat + nbeta2);
+	double p1 = PREDICTOR_INVERSE_LINK_NO_SCALE(p1_intern);
+	double p2 = PREDICTOR_INVERSE_LINK_NO_SCALE(p2_intern);
+	double w1 = dat[BINOMIALMIX_NBETA+2];
+	double w2 = dat[BINOMIALMIX_NBETA+3];
+	double p12 = w1 * p1 + w2 * p2;
+	double w3 = 1.0 - (w1 + w2);
 	assert(w3 >= 0 && w3 <= 1.0);
 
-	double off = OFFSET(idx);
+	double offset = off + beta9 * dat[BINOMIALMIX_NBETA+1];
 	if (m > 0) {
 		if (ISZERO(y)) {
 #pragma omp simd
 			for (int i = 0; i < m; i++) {
-				double p = p12 + w3 * PREDICTOR_INVERSE_LINK(x[i] + off);
+				double p = p12 + w3 * PREDICTOR_INVERSE_LINK(x[i] + offset);
 				logll[i] = normc + ny * LOG_1mp(p);
 			}
 		} else if (ISZERO(ny)) {
 #pragma omp simd
 			for (int i = 0; i < m; i++) {
-				double p = p12 + w3 * PREDICTOR_INVERSE_LINK(x[i] + off);
+				double p = p12 + w3 * PREDICTOR_INVERSE_LINK(x[i] + offset);
 				logll[i] = normc + y * LOG_p(p);
 			}
 		} else {
 #pragma omp simd
 			for (int i = 0; i < m; i++) {
-				double p = p12 + w3 * PREDICTOR_INVERSE_LINK(x[i] + off);
+				double p = p12 + w3 * PREDICTOR_INVERSE_LINK(x[i] + offset);
 				logll[i] = normc + y * LOG_p(p) + ny * LOG_1mp(p);
 			}
 		}
 	} else {
 		double *yy = (y_cdf ? y_cdf : &y);
 		for (int i = 0; i < -m; i++) {
-			double p = p12 + w3 * PREDICTOR_INVERSE_LINK(x[i] + off);
+			double p = p12 + w3 * PREDICTOR_INVERSE_LINK(x[i] + offset);
 			logll[i] = gsl_cdf_binomial_P((unsigned int) *yy, p, (unsigned int) n);
 		}
 	}

inlaprog/src/inla-parse.c

@@ -4883,8 +4883,6 @@ int inla_parse_data(inla_tp *mb, dictionary *ini, int sec)
 	case L_BINOMIALMIX:
 	{
 		const int nbeta = BINOMIALMIX_NBETA;
-		assert(GSL_IS_EVEN(BINOMIALMIX_NBETA));
-
 		for (i = 0; i < nbeta; i++) {
 			GMRFLib_sprintf(&ctmp, "FIXED%1d", i);
 			iniparser_getstring(ini, inla_string_join(secname, ctmp), NULL);

inlaprog/src/inla-testit.c

@@ -2644,7 +2644,7 @@ int testit(int argc, char **argv)
 		P(n);
 		P(m);
 		double sum = 0.0, sum1 = 0.0, sum2 = 0.0, sum0 = 0.0;
-		double tref[3] = {0, 0, 0};
+		double tref[3] = { 0, 0, 0 };
 		double work[n];
 
 		for (int k = 0; k < m; k++) {
@@ -2668,7 +2668,7 @@ int testit(int argc, char **argv)
 		P((sum1 - sum2) / (sum1 + sum2));
 		P((sum - sum2) / (sum + sum2));
 		sum0 = sum;
-		
+
 		tref[0] = tref[1] = tref[2] = sum = sum1 = sum2 = 0.0;
 		for (int k = 0; k < m; k++) {
 			tref[0] -= GMRFLib_timer();
@@ -4469,12 +4469,12 @@ int testit(int argc, char **argv)
 
 #if !defined(INLA_WITH_MKL)
 		mkl = 0;
-#endif		
+#endif
 
 		P(n);
 		P(m);
 		P(mkl);
-		
+
 		for (int i = 0; i < n + 1; i++) {
 			x[i] = GMRFLib_uniform();
 		}

inlaprog/src/inla.c

@@ -121,7 +121,7 @@
 #define RCPOISSON_MAXTHETA (10L)
 #define TPOISSON_MAXTHETA (10L)
 #define OCCUPANCY_MAXTHETA (10L)
-#define BINOMIALMIX_NBETA (8L)
+#define BINOMIALMIX_NBETA (9L)
 #define L_FL_NC (9L)
 
 G_tp G = { 1, INLA_MODE_DEFAULT, 4.0, 0.5, 2, 0, GMRFLib_REORDER_DEFAULT, 0, 0 };
@@ -155,6 +155,7 @@ int R_load_INLA = 0;
 #define OFFSET3(idx_) mb->offset[idx_]
 
 #define LINK_INIT							\
+	double off = OFFSET(idx);					\
 	double *_link_covariates = NULL;				\
 	Link_param_tp *predictor_invlinkfunc_arg = (Link_param_tp *) (ds->predictor_invlinkfunc_arg[idx]); \
 	if (ds->link_covariates) {					\
@@ -166,12 +167,15 @@ int R_load_INLA = 0;
 		_lp_scale = ds->lp_scale_beta[(int)ds->lp_scale[idx]][thread_id][0]; \
 	}
 
+
 #define LINK_END Free(_link_covariates)
 #define PREDICTOR_SCALE _lp_scale
 #define PREDICTOR_LINK_EQ(_fun) (ds->predictor_invlinkfunc == (_fun))
 #define PREDICTOR_SIMPLE_LINK_EQ(_fun) (ds->data_observations.link_simple_invlinkfunc ==  (_fun))
 #define PREDICTOR_INVERSE_LINK(xx_)					\
 	ds->predictor_invlinkfunc(thread_id, _lp_scale * (xx_), MAP_FORWARD, (void *)predictor_invlinkfunc_arg, _link_covariates)
+#define PREDICTOR_INVERSE_LINK_NO_SCALE(xx_)				\
+	ds->predictor_invlinkfunc(thread_id, (xx_), MAP_FORWARD, (void *)predictor_invlinkfunc_arg, _link_covariates)
 #define PREDICTOR_INVERSE_IDENTITY_LINK(xx_) (_lp_scale * (xx_))
 #define PREDICTOR_LINK(xx_)						\
 	(ds->predictor_invlinkfunc(thread_id, (xx_), MAP_BACKWARD, (void *)predictor_invlinkfunc_arg, _link_covariates) / _lp_scale)

inlaprog/src/inla.h

@@ -2285,7 +2285,7 @@ int loglikelihood_logperiodogram(int thread_id, double *logll, double *x, int m,
 int loglikelihood_mgamma(int thread_id, double *logll, double *x, int m, int idx, double *x_vec, double *y_cdf, void *arg, char **arg_str);
 int loglikelihood_mgammasurv(int thread_id, double *logll, double *x, int m, int idx, double *x_vec, double *y_cdf, void *arg, char **arg_str);
 int loglikelihood_mix_core(int thread_id, double *logll, double *x, int m, int idx, double *x_vec, double *y_cdf, void *arg,
-			   int (*quadrature)(int, double **, double **, int *, void *), int(*simpson)(int, double **, double **, int *, void *),
+			   int (*quadrature)(int, double **, double **, int *, void *), int (*simpson)(int, double **, double **, int *, void *),
 			   char **arg_str);
 int loglikelihood_mix_loggamma(int thread_id, double *logll, double *x, int m, int idx, double *x_vec, double *y_cdf, void *arg, char **arg_str);
 int loglikelihood_mix_mloggamma(int thread_id, double *logll, double *x, int m, int idx, double *x_vec, double *y_cdf, void *arg, char **arg_str);

inlaprog/src/my.c

@@ -239,7 +239,7 @@ double my_betabinomial_helper4(int n, double a, double *work)
 #pragma omp simd
 	for (int i = 0; i < nn; i++) {
 		double aa = i * roll + a;
-		work[i] = aa * (aa + 1) * (aa + 2) * (aa + 3); 
+		work[i] = aa * (aa + 1) * (aa + 2) * (aa + 3);
 	}
 
 	GMRFLib_log(nn, work, work);
@@ -329,7 +329,7 @@ double my_betabinomial_helper_core(int n, double a, double *work, int roll)
 		double aa = i * roll + a;
 		double s = 1.0;
 #pragma omp simd reduction(*: s)
-		for(int j = 0; j < roll; j++) {
+		for (int j = 0; j < roll; j++) {
 			s *= (aa + j);
 		}
 		work[i] = s;

r-inla.org/doc/hyper/likelihood/binomialmix.tex

@@ -116,6 +116,20 @@
 	 	 	\item[to.theta] \verb!function(x) x!
 	 	 	\item[from.theta] \verb!function(x) x!
 	 	 \end{description}
+	 	\item[theta9]\ 
+	 	 \begin{description}
+	 	 	\item[hyperid] \verb!56559!
+	 	 	\item[name] \verb!beta9!
+	 	 	\item[short.name] \verb!beta9!
+	 	 	\item[output.name] \verb!beta9 for binomialmix observations!
+	 	 	\item[output.name.intern] \verb!beta9 for binomialmix observations!
+	 	 	\item[initial] \verb!0!
+	 	 	\item[fixed] \verb!FALSE!
+	 	 	\item[prior] \verb!normal!
+	 	 	\item[param] \verb!0 100!
+	 	 	\item[to.theta] \verb!function(x) x!
+	 	 	\item[from.theta] \verb!function(x) x!
+	 	 \end{description}
 	 \end{description}
 	\item[status] \verb!experimental!
 	\item[survival] \verb!FALSE!

r-inla.org/doc/likelihood/binomialmix.tex

@@ -16,29 +16,32 @@ \subsubsection*{Details}
 \begin{displaymath}
     p = w_1 p_1 + w_2 p_2 + w_3 p_3
 \end{displaymath}
-where $w_1 + w_2 + w_3 = 1$ are the positive weights, and
+where $w_1 + w_2 + w_3 = 1$ are the positive weights,
 \begin{displaymath}
-    \text{logit}(p_3) = \eta
+    \text{logit}(p_1) = \beta_{1} z_{1} + \beta_{2} z_{2} +
+    \beta_{3} z_{3} + \beta_4 z_4 + \beta_9 z_9,
 \end{displaymath}
-where the linear predictor $\eta$ is defined in the formula. Further,
 \begin{displaymath}
-    \text{logit}(p_1) = \beta_{1} z_{1} + \beta_{2} z_{2} +
-    \beta_{3} z_{3} + \beta_4 z_4
+    \text{logit}(p_2) = \beta_{5} z_{5} +
+    \beta_{6} z_{6} + \beta_7 z_7 + \beta_8 z_8 + \beta_9 z_{10}
 \end{displaymath}
 and
 \begin{displaymath}
-    \text{logit}(p_2) = \beta_{5} z_{5} +
-    \beta_{6} z_{6} + \beta_7 z_7 + \beta_8 z_8
+    \text{logit}(p_3) = \eta + \beta_9 z_{11}
 \end{displaymath}
-for fixed covariates $z_{i}$. 
+for fixed covariates $z_{i}$.
+The linear predictor $\eta$ is defined in the formula.
+
+\textbf{Note:} $\beta_9$ is the \emph{same} variable in the three
+expressions.
 
 \subsection*{Link-function}
 
 The link-function is given as usual, and they are all equal.
 
 \subsection*{Hyperparameters}
 
-The eight regression coefficients $\beta_{i}$ are treated as
+The nine regression coefficients $\beta_{i}$ are treated as
 hyperparameters. 
 
 \subsection*{Specification}
@@ -47,7 +50,7 @@ \subsection*{Specification}
 \item \texttt{family="binomialmix"}
 \item Required arguments: A $n\times 2$ matrix $Y$ with the
 observations and the number of trials $s$,  $Y=(y, s)$,
-    a $n\times 8$ matrix $Z$ with the covariates $Z=c(z_1, \ldots, z_8)$,
+    a $n\times 11$ matrix $Z$ with the covariates $Z=c(z_1, \ldots, z_{11})$,
     and a $n\times 2$ matrix $W$ with weights $W=(w_1, w_2)$. The
     \texttt{inla.mdata} is used as
 \begin{verbatim}

r-inla.org/doc/likelihood/example-binomialmix.R

@@ -1,37 +1,41 @@
-n <- 30000
+n <- 10^5
+## high value of 'size' makes it easy to check the likelihood implementation
+size <- sample(100:200, n, replace = TRUE)
 
-## this makes it to easy to its just to check the likelihood
-## implementation
-size <- sample(10:20, n, replace = TRUE)
+beta.p1 <- rnorm(4, sd = 0.5)
+beta.p2 <- rnorm(4, sd = 0.5)
+beta9 <- rnorm(1, sd = 0.5)
+beta <- c(beta.p1, beta.p2, beta9)
 
-beta1 <- rnorm(4, sd = 0.2)
-beta2 <- rnorm(4, sd = 0.2)
-beta <- c(beta1, beta2)
-Z <- matrix(NA, n, 8)
+Z <- matrix(NA, n, 11)
 W <- matrix(NA, n, 2)
 Y <- matrix(NA, n, 2)
 
-x <- rnorm(n, sd = 0.1)
-xx <- rnorm(n, sd = 0.4)
-eta <-  1 + x + xx
+x <- rnorm(n, sd = 0.5)
+xx <- rnorm(n, sd = 0.5)
+eta <- numeric(n)
 
 for (i in 1:n) {
-    Z[i, ] <- rnorm(8)
+    Z[i, ] <- rnorm(11)
     w <- c(rbeta(2, 1, 10), rbeta(1, 10, 1))
     w <- w/sum(w)
     W[i, ] <- w[1:2]
 
-    p1 <- inla.link.invlogit(sum(beta1 * Z[i, 1:4]))
-    p2 <- inla.link.invlogit(sum(beta2 * Z[i, 4 + 1:4]))
+    p1 <- inla.link.invlogit(sum(beta.p1 * Z[i, 1:4]) + beta9 * Z[i, 9])
+    p2 <- inla.link.invlogit(sum(beta.p2 * Z[i, 4 + 1:4]) + beta9 * Z[i, 10])
+
+    eta[i] <-  1 + x[i] + xx[i] + beta9 * Z[i, 11]
     p3 <- inla.link.invlogit(eta[i])
+
     p <- w[1] * p1 + w[2] * p2 + w[3] * p3
     Y[i, ] <- c(rbinom(1, size = size[i], prob = p), size[i])
 }
 
 r <- inla(inla.mdata(Y, Z, W) ~ 1 + x + xx, 
           family = "binomialmix",
           data = list(Y = Y, Z = Z, W = W, x = x, xx = xx),
-          verbose = TRUE)
+          verbose = TRUE,
+          control.inla = list(int.strategy = "eb"))
 
 print(round(dig = 4, cbind(estimate = r$summary.fixed[,"mean"], true = 1)))
 print(round(dig = 4, cbind(estimate =  r$summary.hyperpar[,"mean"], true = beta)))