#include #include #include #define PI 3.14159265 //do product recursively //fixed x, n, g, p, but different mu, sigma for different p; need an index vector id[] (i_j in equation 53) for this //remember both k and j start from 1 to n and p, respectively; need to subtract one for index consistency double prod_of(long *n, long *g, long *p, long *k, long *j, long *id, double *x, double *mu, double *sigma){ int i; double prod; long jnext; if ((*j)==1){ prod=exp(0-0.5*(x[((*k)-1)*(*p)+(*j)-1]-mu[((*j)-1)*(*g)+id[(*j)-1]])*(x[((*k)-1)*(*p)+(*j)-1]-mu[((*j)-1)*(*g)+id[(*j)-1]])/sigma[((*j)-1)*(*g)+id[(*j)-1]])/sqrt(2*PI*sigma[((*j)-1)*(*g)+id[(*j)-1]]); } else { jnext = (*j)-1; prod=prod_of(n,g,p,k,&jnext,id,x,mu,sigma)*exp(0-0.5*(x[((*k)-1)*(*p)+(*j)-1]-mu[((*j)-1)*(*g)+id[(*j)-1]])*(x[((*k)-1)*(*p)+(*j)-1]-mu[((*j)-1)*(*g)+id[(*j)-1]])/sigma[((*j)-1)*(*g)+id[(*j)-1]])/sqrt(2*PI*sigma[((*j)-1)*(*g)+id[(*j)-1]]); } return(prod); /*printf("%lf \n", prod);*/ } //make an index vector for each of combination of dimensional components //remember j start from 1 to p; need to subtract one for index consistency long index_of(long *g, long *j, long *id){ long idx; long jnext; if ((*j)==1){ idx=id[(*j)-1]; } else { jnext=(*j)-1; idx=index_of(g,&jnext,id)*(*g)+id[(*j)-1]; } return(idx); } //calculate partial density for diagonal elements double fpcdpart1_k(long *n, long *g, long *p, long *k, double *pi, double *alpha, double *x, double *mu, double *sigma){ double total; int l,jj; long id[(*p)]; total=0; for (l=0; l<(*g); l++){ for (jj=0; jj<(*p); jj++){ id[jj]=l; } total=total+(pi[l]-(*alpha))*prod_of(n, g, p, k, p, id, x, mu, sigma); } return(total); } //calculate partial density for non-diagonal elements //do multiple sums using for and also recursion //remember j start from 1 to p; need to subtract one for index consistency double nest_of_fpcdpart2_k(long *n, long *g, long *p, long *k, long *j, long *id, double *x, double *mu, double *sigma){ double total; int l; long jnext; l=0-1; total=0; if ((*j)<=(*p)){ for(l=0;l<(*g);l++){ id[(*j)-1]=l; jnext=(*j)+1; total=total+nest_of_fpcdpart2_k(n, g, p, k, &jnext, id, x, mu, sigma); } } else { total=prod_of(n, g, p, k, p, id, x, mu, sigma); } /*printf("%d %d %d %lf \n", *k, *j, l, total);*/ return(total); } //calculate uhat which includes alpha void nest_ez_k(long *n, long *g, long *p, long *k, long *j, long *id, double *uhat, double *pi, double *alpha, double *x, double *mu, double *sigma){ double temp1, temp2, temp4, sumpi; int l; long jnext, idx, ii,jj,temp, temp3, id1[(*p)]; fflush(stdin); l=0-1; jj=1; temp=pow((*g),(*p)); if ((*j)<=(*p)) { for (l=0; l<(*g); l++){ id[(*j)-1]=l; jnext=(*j)+1; nest_ez_k(n,g,p,k,&jnext, id, uhat, pi, alpha, x, mu, sigma); } /*for(ii=0; ii<(*p); ii++){ printf("%d ", id[ii]); }*/ } else { temp1=fpcdpart1_k(n, g, p, k, pi, alpha, x, mu, sigma); temp2=nest_of_fpcdpart2_k(n, g, p, k, &jj, id1, x, mu, sigma); temp4=prod_of(n, g, p, k, p, id, x, mu, sigma); /*printf("%lf \n",temp1);*/ temp3=index_of(g,p,id); /*printf("%d", temp3); printf("\n");*/ uhat[((*k)-1)*temp+temp3]=((*alpha)*temp4)/(temp1+(*alpha)*temp2); /*for(ii=0; ii<(*p); ii++){ printf("%d ", id[ii]); } printf("\n");*/ } } double nest_of_mu_jj_num_k(long *n, long *g, long *p, long *k, long *j, long *jj, long *ijj, long *id ,double *uhat, double *pi, double *x, double *mu, double *sigma){ double total; int l; long jnext, temp1, temp; fflush(stdin); temp=pow((*g),(*p)); l=0-1; total=0; if ((*j)<=(*p)){ if ((*j)==(*jj)){ for(l=0;l<(*g);l++){ if (l==(*ijj)){ id[(*j)-1]=l; jnext=(*j)+1; total=total+nest_of_mu_jj_num_k(n, g, p, k, &jnext, jj, ijj, id, uhat, pi, x, mu, sigma); /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/ } else{ total=total+0; /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/ } } } else{ for(l=0;l<(*g);l++){ id[(*j)-1]=l; jnext=(*j)+1; total=total+nest_of_mu_jj_num_k(n, g, p, k, &jnext, jj, ijj, id, uhat, pi, x, mu, sigma); /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/ } } } else { if ((*j)!=(*jj)){ temp1=index_of(g,p,id); total=uhat[((*k)-1)*temp+temp1]*x[((*k)-1)*(*p)+(*jj)-1]; } /*else { if (id[(*p)-1]==(*ijj)) { temp1=index_of(g,p,id); total=uhat[((*k)-1)*temp+temp1]*x[((*k)-1)*(*p)+(*jj)-1]; } }*/ } /*printf("%d %d %d %lf \n", *k, *j, l, total);*/ return(total); } double nest_of_mu_jj_den_k(long *n, long *g, long *p, long *k, long *j, long *jj, long *ijj, long *id, double *uhat, double *pi, double *x, double *mu, double *sigma){ double total; int l; long jnext, temp1, temp; fflush(stdin); temp=pow((*g),(*p)); l=0-1; total=0; if ((*j)<=(*p)){ if ((*j)==(*jj)){ for(l=0;l<(*g);l++){ if (l==(*ijj)){ id[(*j)-1]=l; jnext=(*j)+1; total=total+nest_of_mu_jj_den_k(n, g, p, k, &jnext, jj, ijj, id, uhat, pi, x, mu, sigma); /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/ } else{ total=total+0; /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/ } } } else{ for(l=0;l<(*g);l++){ id[(*j)-1]=l; jnext=(*j)+1; total=total+nest_of_mu_jj_den_k(n, g, p, k, &jnext, jj, ijj, id, uhat, pi, x, mu, sigma); /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/ } } } else { if ((*j)!=(*jj)){ temp1=index_of(g,p,id); total=uhat[((*k)-1)*temp+temp1]; } /*else { if (id[(*p)-1]==(*ijj)) { temp1=index_of(g,p,id); total=uhat[((*k)-1)*temp+temp1]; } }*/ } /*printf("%d %d %d %lf \n", *k, *j, l, total);*/ return(total); } double nest_of_sigma_jj_num_k(long *n, long *g, long *p, long *k, long *j, long *jj, long *ijj, long *id, double *uhat, double *pi, double *x, double *mu, double *sigma){ double total; int l; long jnext, temp1, temp; fflush(stdin); temp=pow((*g),(*p)); l=0-1; total=0; if ((*j)<=(*p)){ if ((*j)==(*jj)){ for(l=0;l<(*g);l++){ if (l==(*ijj)){ id[(*j)-1]=l; jnext=(*j)+1; total=total+nest_of_sigma_jj_num_k(n, g, p, k, &jnext, jj, ijj,id, uhat, pi, x, mu, sigma); /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/ } else{ total=total+0; /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/ } } } else{ for(l=0;l<(*g);l++){ id[(*j)-1]=l; jnext=(*j)+1; total=total+nest_of_sigma_jj_num_k(n, g, p, k, &jnext, jj, ijj,id, uhat, pi, x, mu, sigma); /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/ } } } else { if ((*j)!=(*jj)){ temp1=index_of(g,p,id); total=uhat[((*k)-1)*temp+temp1]*(x[((*k)-1)*(*p)+(*jj)-1]-mu[((*jj)-1)*(*g)+(*ijj)])*(x[((*k)-1)*(*p)+(*jj)-1]-mu[((*jj)-1)*(*g)+(*ijj)]); } /*else { if (id[(*p)-1]==(*ijj)) { temp1=index_of(g,p,id); total=uhat[((*k)-1)*temp+temp1]*(x[((*k)-1)*(*p)+(*jj)-1]-mu[((*jj)-1)*(*g)+(*ijj)])*(x[((*k)-1)*(*p)+(*jj)-1]-mu[((*jj)-1)*(*g)+(*ijj)]); } }*/ } /*printf("%d %d %d %lf \n", *k, *j, l, total);*/ return(total); } void em_pdim_pcd_exchange(long *n, long *g, long *p, long *j, long *id, double *uhat, double *pi, double *alpha, double *x, double *mu, double *sigma, double *loglik,double *tol, long *restriction, long *constrain, long *iteration, long *convergence){ long h,l,jj,ij,jjj,ijj,idx,*id1,*id2,power; int i, ii,flag, iter; double *Lg,*logLg,loglik2,sumpi,temppi_den, tempmu_num_sum, tempmu_den_sum, tempsigma_num_sum, *tempmu_num_jj,*tempmu_den_jj,*tempsigma_num_jj, *fpcd; power=pow((*g),(*p)); id1 = (long *) malloc((*p) * sizeof(long)); id2 = (long *) malloc((*p) * sizeof(long)); tempmu_num_jj= (double *) malloc((*n) * sizeof(double)); tempmu_den_jj = (double *) malloc((*n)* sizeof(double)); tempsigma_num_jj = (double *) malloc((*n) * sizeof(double)); fpcd = (double *) malloc((*n) * sizeof(double)); logLg = (double *) malloc((*n) * sizeof(double)); /*theoretical restriction on the null*/ if((*restriction)>0){ /*initialization*/ flag=1; (*loglik) = 0; for(h=1;h<=(*n);h++){ for(i=0;i<(*p);i++){ id[i]=0; } fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi, alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma); logLg[h-1] = log(fpcd[h-1]); (*loglik)= (*loglik) + logLg[h-1]; // printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]); } printf("%lf \n", *loglik); /*iteration*/ iter = 0; while(flag>0 & iter<(*iteration)){ printf("%d \n",iter); /*E step for z*/ for(h=1;h<=(*n);h++){ for(i=0;i<(*p);i++){ id[i]=0; } nest_ez_k(n, g, p, &h, j, id, uhat, pi,alpha, x, mu, sigma); } for (l=0;l0.0){ mu[jj*(*g)+ij]=0.0; } } if(constrain[jj*(*g)+ij]==0){ mu[jj*(*g)+ij] = 0.0; sigma[jj*(*g)+ij] = 1.0; } if(constrain[jj*(*g)+ij]==0+1){ if(mu[jj*(*g)+ij]<0.0){ mu[jj*(*g)+ij]=0.0; } } } } /*check convergence*/ loglik2 = (*loglik); (*loglik) = 0; for(h=1;h<=(*n);h++){ for(i=0;i<(*p);i++){ id[i]=0; } fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi,alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma); logLg[h-1] = log(fpcd[h-1]); (*loglik)= (*loglik) + logLg[h-1]; // printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]); } if(fabs((*loglik)-loglik2)<(*tol)){ flag = 0; } iter = iter + 1; } (*convergence) = 0; if(iter<(*iteration)){ (*convergence) = 1; } for (l=0;l<(*g);l++){ printf("%d %lf \n", l, pi[l]); } printf("%lf \n", *alpha); for (l=0;l<(*g)*(*p);l++){ printf("%d %lf %lf \n", l, mu[l],sigma[l]); } } /*no restriction on the null*/ else{ /*initialization*/ flag=1; (*loglik) = 0; for(h=1;h<=(*n);h++){ for(i=0;i<(*p);i++){ id[i]=0; } fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi,alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma); logLg[h-1] = log(fpcd[h-1]); (*loglik)= (*loglik) + logLg[h-1]; printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]); } /*printf("%lf \n", *loglik);/ /*iteration*/ iter = 0; while(flag>0 & iter<(*iteration)){ /*E step for z*/ for(h=1;h<=(*n);h++){ for(i=0;i<(*p);i++){ id[i]=0; } nest_ez_k(n, g, p, &h, j, id, uhat, pi, alpha, x, mu, sigma); } for (l=0;l0){ /*initialization*/ flag=1; (*loglik) = 0; for(h=1;h<=(*n);h++){ for(i=0;i<(*p);i++){ id[i]=0; } fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi, alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma); logLg[h-1] = log(fpcd[h-1]); (*loglik)= (*loglik) + logLg[h-1]; // printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]); } printf("%lf \n", *loglik); /*iteration*/ iter = 0; while(flag>0 & iter<(*iteration)){ printf("%d \n",iter); /*M step for pi*/ temppi_den=0; for(l=0;l0.0){ mu[jj*(*g)+ij]=0.0; } } if(constrain[jj*(*g)+ij]==0){ mu[jj*(*g)+ij] = 0.0; sigma[jj*(*g)+ij] = 1.0; } if(constrain[jj*(*g)+ij]==0+1){ if(mu[jj*(*g)+ij]<0.0){ mu[jj*(*g)+ij]=0.0; } } } } /*check convergence*/ loglik2 = (*loglik); (*loglik) = 0; for(h=1;h<=(*n);h++){ for(i=0;i<(*p);i++){ id[i]=0; } fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi,alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma); logLg[h-1] = log(fpcd[h-1]); (*loglik)= (*loglik) + logLg[h-1]; // printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]); } if(fabs((*loglik)-loglik2)<(*tol)){ flag = 0; } iter = iter + 1; } (*convergence) = 0; if(iter<(*iteration)){ (*convergence) = 1; } for (l=0;l<(*g);l++){ printf("%d %lf \n", l, pi[l]); } printf("%lf \n", *alpha); for (l=0;l<(*g)*(*p);l++){ printf("%d %lf %lf \n", l, mu[l],sigma[l]); } } /*no restriction on the null*/ else{ /*initialization*/ flag=1; (*loglik) = 0; for(h=1;h<=(*n);h++){ for(i=0;i<(*p);i++){ id[i]=0; } fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi,alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma); logLg[h-1] = log(fpcd[h-1]); (*loglik)= (*loglik) + logLg[h-1]; printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]); } /*printf("%lf \n", *loglik);/ /*iteration*/ iter = 0; while(flag>0 & iter<(*iteration)){ /*M step for pi*/ temppi_den=0; for(l=0;l