- for (int r=0; r<k; r++)
- gam[mi(i,r,n,k)] /= sumGamI;
+ logGam[r] = logGam[r] - maxLogGam; //adjust without changing proportions
+ gam[mi(i,r,n,k)] = exp(logGam[r]); //gam[i, ] <- exp(logGam)
+ norm_fact += gam[mi(i,r,n,k)]; //norm_fact <- sum(gam[i, ])