/*=======================================================================================*/ // Author: Ryan Martin // C subroutine for computing PR estimates and marginal likelihood for finite mixtures. // See corresponding paper for details of the method. // These subroutines are to be called by "sasa.R" -- an R front-end function. /*=======================================================================================*/ #include #include #include #include #include double vec_sum(double *x, int len) { int i; double s = 0.0; for(i = 0; i < len; i++) s += x[i]; return s; } void PR_finite(double *lik, int *nobs, int *ngrid, double *w, double *f, double *logmix, int *nperm, int *perm) { int i, j, k; double sum_post, L = 0.0, LL; double *post = (double *)R_alloc(ngrid[0], sizeof(double)); double *F = (double *)R_alloc(ngrid[0], sizeof(double)); double *ff = (double *)R_alloc(ngrid[0], sizeof(double)); for(k = 0; k < nperm[0]; k++) { LL = logmix[0]; for(j = 0; j < ngrid[0]; j++) ff[j] = f[j]; for(i = 0; i < nobs[0]; i++) { for(j = 0; j < ngrid[0]; j++) post[j] = ff[j] * lik[perm[i + nobs[0] * k] + nobs[0] * j]; sum_post = vec_sum(post, ngrid[0]); LL += log(sum_post); for(j = 0; j < ngrid[0]; j++) ff[j] += w[i] * (post[j] / sum_post - ff[j]); } L += (LL - L) / (k + 1); for(j = 0; j < ngrid[0]; j++) F[j] += (ff[j] - F[j]) / (k + 1); } logmix[0] = -L; for(j = 0; j < ngrid[0]; j++) f[j] = F[j]; }