#include <math.h>

void compute_zstats(float **zmean,
		    float **zvar_inv,
		    float *zcorprod,
		    float **wta,
		    float **wtb,
		    float **f,
		    float **xmean,
		    float **xvar,
		    float *xprior,
		    float *n, float *q, float *n_var, int nmodes, int ndim)
{
	int i, j;

	for (i = 0; i < nmodes; ++i) {
		zcorprod[i] = log(xprior[i]);
		for (j = 0; j < ndim; ++j) {
			zmean[i][j] =
			    wta[i][j] * (xmean[i][j] + q[j]) +
			    wtb[i][j] * n[j] + f[i][j];
			zvar_inv[i][j] =
			    1.0 / (wta[i][j] * wta[i][j] * xvar[i][j] +
				   wtb[i][j] * wtb[i][j] * n_var[j]);
			zcorprod[i] += 0.5 * log(zvar_inv[i][j]);
		}
	}
	return;
}
