{
gsl_matrix* U = gsl_matrix_alloc(dim,dim);
gsl_matrix* V = gsl_matrix_alloc(dim,dim);
{
gsl_matrix* U = gsl_matrix_alloc(dim,dim);
gsl_matrix* V = gsl_matrix_alloc(dim,dim);
//U,S,V = SVD of matrix
gsl_linalg_SV_decomp(U, V, S, work);
gsl_vector_free(work);
//U,S,V = SVD of matrix
gsl_linalg_SV_decomp(U, V, S, work);
gsl_vector_free(work);
dotProduct += V->data[i*dim+j] * (S->data[j] > EPS ? 1.0/S->data[j] : 0.0) * U->data[ii*dim+j];
inverse[i*dim+ii] = dotProduct;
}
dotProduct += V->data[i*dim+j] * (S->data[j] > EPS ? 1.0/S->data[j] : 0.0) * U->data[ii*dim+j];
inverse[i*dim+ii] = dotProduct;
}
- const Real* Pi, // parametre de proportion
- const Real* Rho, // parametre initial de variance renormalisé
- Int mini, // nombre minimal d'itérations dans l'algorithme EM
- Int maxi, // nombre maximal d'itérations dans l'algorithme EM
- const Real* X, // régresseurs
- const Real* Y, // réponse
- Real tau, // seuil pour accepter la convergence
- const Int* rank, // vecteur des rangs possibles
+ const Real* Pi, // parametre de proportion
+ const Real* Rho, // parametre initial de variance renormalisé
+ int mini, // nombre minimal d'itérations dans l'algorithme EM
+ int maxi, // nombre maximal d'itérations dans l'algorithme EM
+ const Real* X, // régresseurs
+ const Real* Y, // réponse
+ Real tau, // seuil pour accepter la convergence
+ const int* rank, // vecteur des rangs possibles
- Real* phi, // parametre de moyenne renormalisé, calculé par l'EM
- Real* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres
+ Real* phi, // parametre de moyenne renormalisé, calculé par l'EM
+ Real* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres
- mwSize n, // taille de l'echantillon
- mwSize p, // nombre de covariables
- mwSize m, // taille de Y (multivarié)
- mwSize k) // nombre de composantes
+ int n, // taille de l'echantillon
+ int p, // nombre de covariables
+ int m, // taille de Y (multivarié)
+ int k) // nombre de composantes
{
// Allocations, initializations
Real* Phi = (Real*)calloc(p*m*k,sizeof(Real));
{
// Allocations, initializations
Real* Phi = (Real*)calloc(p*m*k,sizeof(Real));
Real* Yr = (Real*)malloc(n*m*sizeof(Real));
Real* tXrXr = (Real*)malloc(p*p*sizeof(Real));
Real* tXrYr = (Real*)malloc(p*m*sizeof(Real));
Real* Yr = (Real*)malloc(n*m*sizeof(Real));
Real* tXrXr = (Real*)malloc(p*p*sizeof(Real));
Real* tXrYr = (Real*)malloc(p*m*sizeof(Real));
gsl_matrix* matrixE = gsl_matrix_alloc(m, m);
gsl_permutation* permutation = gsl_permutation_alloc(m);
gsl_matrix* V = gsl_matrix_alloc(m,m);
gsl_matrix* matrixE = gsl_matrix_alloc(m, m);
gsl_permutation* permutation = gsl_permutation_alloc(m);
gsl_matrix* V = gsl_matrix_alloc(m,m);
gsl_vector* work = gsl_vector_alloc(m);
//Initialize class memberships (all elements in class 0; TODO: randomize ?)
gsl_vector* work = gsl_vector_alloc(m);
//Initialize class memberships (all elements in class 0; TODO: randomize ?)
- for (mwSize j=0; j<p; j++)
- Xr[cardClustR*p+j] = X[i*p+j];
- for (mwSize j=0; j<m; j++)
- Yr[cardClustR*m+j] = Y[i*m+j];
+ for (int j=0; j<p; j++)
+ Xr[mi(cardClustR,j,n,p)] = X[mi(i,j,n,p)];
+ for (int j=0; j<m; j++)
+ Yr[mi(cardClustR,j,n,m)] = Y[mi(i,j,n,m)];
- for (mwSize u=0; u<cardClustR; u++)
- dotProduct += Xr[u*p+j] * Xr[u*p+jj];
- tXrXr[j*p+jj] = dotProduct;
+ for (int u=0; u<cardClustR; u++)
+ dotProduct += Xr[mi(u,j,n,p)] * Xr[mi(u,jj,n,p)];
+ tXrXr[mi(j,jj,p,p)] = dotProduct;
Real* invTXrXr = pinv(tXrXr, p);
// Compute tXrYr = t(Xr) * Yr
Real* invTXrXr = pinv(tXrXr, p);
// Compute tXrYr = t(Xr) * Yr
- for (mwSize u=0; u<cardClustR; u++)
- dotProduct += Xr[u*p+j] * Yr[u*m+jj];
- tXrYr[j*m+jj] = dotProduct;
+ for (int u=0; u<cardClustR; u++)
+ dotProduct += Xr[mi(u,j,n,p)] * Yr[mi(u,jj,n,m)];
+ tXrYr[mi(j,jj,p,m)] = dotProduct;
- for (mwSize u=0; u<p; u++)
- dotProduct += invTXrXr[j*p+u] * tXrYr[u*m+jj];
- matrixM->data[j*m+jj] = dotProduct;
+ for (int u=0; u<p; u++)
+ dotProduct += invTXrXr[mi(j,u,p,p)] * tXrYr[mi(u,jj,p,m)];
+ matrixM->data[j*m+jj] = dotProduct;
//U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
gsl_linalg_SV_decomp(matrixM, V, S, work);
//U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
gsl_linalg_SV_decomp(matrixM, V, S, work);
//Compute phi(:,:,r) = hatBetaR * Rho(:,:,r)
//Compute phi(:,:,r) = hatBetaR * Rho(:,:,r)
- for (mwSize u=0; u<m; u++)
- dotProduct += hatBetaR[j*m+u] * Rho[u*m*k+jj*k+r];
- phi[j*m*k+jj*k+r] = dotProduct;
+ for (int u=0; u<m; u++)
+ dotProduct += hatBetaR[mi(j,u,p,m)] * Rho[ai(u,jj,r,m,m,k)];
+ phi[ai(j,jj,r,p,m,k)] = dotProduct;
{
//Compute
//Gam(i,r) = Pi(r) * det(Rho(:,:,r)) * exp( -1/2 * (Y(i,:)*Rho(:,:,r) - X(i,:)...
{
//Compute
//Gam(i,r) = Pi(r) * det(Rho(:,:,r)) * exp( -1/2 * (Y(i,:)*Rho(:,:,r) - X(i,:)...
- // *phi(:,:,r)) * transpose( Y(i,:)*Rho(:,:,r) - X(i,:)*phi(:,:,r) ) );
+ //*phi(:,:,r)) * transpose( Y(i,:)*Rho(:,:,r) - X(i,:)*phi(:,:,r) ) );
//split in several sub-steps
//compute det(Rho(:,:,r)) [TODO: avoid re-computations]
//split in several sub-steps
//compute det(Rho(:,:,r)) [TODO: avoid re-computations]
- for (mwSize jj=0; jj<m; jj++)
- matrixE->data[j*m+jj] = Rho[j*m*k+jj*k+r];
+ for (int jj=0; jj<m; jj++)
+ matrixE->data[j*m+jj] = Rho[ai(j,jj,r,m,m,k)];
}
gsl_linalg_LU_decomp(matrixE, permutation, &signum);
Real detRhoR = gsl_linalg_LU_det(matrixE, signum);
}
gsl_linalg_LU_decomp(matrixE, permutation, &signum);
Real detRhoR = gsl_linalg_LU_det(matrixE, signum);
//compute Y(i,:)*Rho(:,:,r)
//compute Y(i,:)*Rho(:,:,r)
- for (mwSize u=0; u<m; u++)
- YiRhoR[j] += Y[i*m+u] * Rho[u*m*k+j*k+r];
+ for (int u=0; u<m; u++)
+ YiRhoR[j] += Y[mi(i,u,n,m)] * Rho[ai(u,j,r,m,m,k)];
//compute X(i,:)*phi(:,:,r)
//compute X(i,:)*phi(:,:,r)
- for (mwSize u=0; u<p; u++)
- XiPhiR[j] += X[i*p+u] * phi[u*m*k+j*k+r];
+ for (int u=0; u<p; u++)
+ XiPhiR[j] += X[mi(i,u,n,p)] * phi[ai(u,j,r,p,m,k)];
//compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
Real dotProduct = 0.0;
//compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
Real dotProduct = 0.0;
dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
Real logGamIR = log(Pi[r]) + log(detRhoR) - 0.5*dotProduct;
dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
Real logGamIR = log(Pi[r]) + log(detRhoR) - 0.5*dotProduct;
//Z(i) = index of max (gam(i,:))
if (logGamIR > maxLogGamIR)
{
//Z(i) = index of max (gam(i,:))
if (logGamIR > maxLogGamIR)
{
// Assign output variable LLF
*LLF = -invN * sumLogLLF2;
//newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi))));
Real newDeltaPhi = 0.0;
// Assign output variable LLF
*LLF = -invN * sumLogLLF2;
//newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi))));
Real newDeltaPhi = 0.0;
- Real tmpDist = fabs(phi[j*m*k+jj*k+r]-Phi[j*m*k+jj*k+r])
- / (1.0+fabs(phi[j*m*k+jj*k+r]));
+ Real tmpDist = fabs(phi[ai(j,jj,r,p,m,k)]-Phi[ai(j,jj,r,p,m,k)])
+ / (1.0+fabs(phi[ai(j,jj,r,p,m,k)]));
- for (mwSize r=0; r<k; r++)
- Phi[j*m*k+jj*k+r] = phi[j*m*k+jj*k+r];
+ for (int r=0; r<k; r++)
+ Phi[ai(jj,j,r,p,m,k)] = phi[ai(jj,j,r,p,m,k)];