diff options
Diffstat (limited to 'src/logistic.cpp')
| -rw-r--r-- | src/logistic.cpp | 746 |
1 files changed, 367 insertions, 379 deletions
diff --git a/src/logistic.cpp b/src/logistic.cpp index f9edc68..2dd0402 100644 --- a/src/logistic.cpp +++ b/src/logistic.cpp @@ -1,15 +1,16 @@ -#include <stdio.h> -#include <math.h> +#include <gsl/gsl_linalg.h> #include <gsl/gsl_matrix.h> -#include <gsl/gsl_rng.h> #include <gsl/gsl_multimin.h> +#include <gsl/gsl_rng.h> #include <gsl/gsl_sf.h> -#include <gsl/gsl_linalg.h> +#include <math.h> +#include <stdio.h> + #include "logistic.h" // I need to bundle all the data that goes to the function to optimze // together. -typedef struct{ +typedef struct { gsl_matrix_int *X; gsl_vector_int *nlev; gsl_vector *y; @@ -18,13 +19,9 @@ typedef struct{ double lambdaL2; } fix_parm_mixed_T; -double fLogit_mixed(gsl_vector *beta, - gsl_matrix_int *X, - gsl_vector_int *nlev, - gsl_matrix *Xc, - gsl_vector *y, - double lambdaL1, - double lambdaL2) { +double fLogit_mixed(gsl_vector *beta, gsl_matrix_int *X, gsl_vector_int *nlev, + gsl_matrix *Xc, gsl_vector *y, double lambdaL1, + double lambdaL2) { int n = y->size; int npar = beta->size; double total = 0; @@ -33,57 +30,56 @@ double fLogit_mixed(gsl_vector *beta, // Changed loop start at 1 instead of 0 to avoid regularization of // beta_0*\/ // #pragma omp parallel for reduction (+:total) - for(int i = 1; i < npar; ++i) - total += beta->data[i]*beta->data[i]; - total = (-total*lambdaL2/2); + for (int i = 1; i < npar; ++i) + total += beta->data[i] * beta->data[i]; + total = (-total * lambdaL2 / 2); // #pragma omp parallel for reduction (+:aux) - for(int i = 1; i < npar; ++i) - aux += (beta->data[i]>0 ? beta->data[i] : -beta->data[i]); - total = total-aux*lambdaL1; + for (int i = 1; i < npar; ++i) + aux += (beta->data[i] > 0 ? beta->data[i] : -beta->data[i]); + total = total - aux * lambdaL1; // #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y) // #reduction (+:total) - for(int i = 0; i < n; ++i) { - double Xbetai=beta->data[0]; - int iParm=1; - for(int k = 0; k < X->size2; ++k) { - if(gsl_matrix_int_get(X,i,k)>0) - Xbetai+=beta->data[gsl_matrix_int_get(X,i,k)-1+iParm]; - iParm+=nlev->data[k]-1; + for (int i = 0; i < n; ++i) { + double Xbetai = beta->data[0]; + int iParm = 1; + for (int k = 0; k < X->size2; ++k) { + if (gsl_matrix_int_get(X, i, k) > 0) + Xbetai += beta->data[gsl_matrix_int_get(X, i, k) - 1 + iParm]; + iParm += nlev->data[k] - 1; } - for(int k = 0; k < (Xc->size2); ++k) - Xbetai+= gsl_matrix_get(Xc,i,k)*beta->data[iParm++]; - total += y->data[i]*Xbetai-gsl_sf_log_1plusx(gsl_sf_exp(Xbetai)); + for (int k = 0; k < (Xc->size2); ++k) + Xbetai += gsl_matrix_get(Xc, i, k) * beta->data[iParm++]; + total += y->data[i] * Xbetai - gsl_sf_log_1plusx(gsl_sf_exp(Xbetai)); } return -total; } void logistic_mixed_pred(gsl_vector *beta, // Vector of parameters - // length = 1 + Sum_k(C_k -1) - gsl_matrix_int *X, // Matrix Nobs x K - gsl_vector_int *nlev, // Vector with number categories - gsl_matrix *Xc, // Continuous covariates matrix: - // obs x Kc (NULL if not used). - gsl_vector *yhat){ // Vector of prob. predicted by - // the logistic - for(int i = 0; i < X->size1; ++i) { - double Xbetai=beta->data[0]; - int iParm=1; - for(int k = 0; k < X->size2; ++k) { - if(gsl_matrix_int_get(X,i,k)>0) - Xbetai+=beta->data[gsl_matrix_int_get(X,i,k)-1+iParm]; - iParm+=nlev->data[k]-1; + // length = 1 + Sum_k(C_k -1) + gsl_matrix_int *X, // Matrix Nobs x K + gsl_vector_int *nlev, // Vector with number categories + gsl_matrix *Xc, // Continuous covariates matrix: + // obs x Kc (NULL if not used). + gsl_vector *yhat) { // Vector of prob. predicted by + // the logistic + for (int i = 0; i < X->size1; ++i) { + double Xbetai = beta->data[0]; + int iParm = 1; + for (int k = 0; k < X->size2; ++k) { + if (gsl_matrix_int_get(X, i, k) > 0) + Xbetai += beta->data[gsl_matrix_int_get(X, i, k) - 1 + iParm]; + iParm += nlev->data[k] - 1; } // Adding the continuous. - for(int k = 0; k < (Xc->size2); ++k) - Xbetai+= gsl_matrix_get(Xc,i,k)*beta->data[iParm++]; - yhat->data[i]=1/(1 + gsl_sf_exp(-Xbetai)); + for (int k = 0; k < (Xc->size2); ++k) + Xbetai += gsl_matrix_get(Xc, i, k) * beta->data[iParm++]; + yhat->data[i] = 1 / (1 + gsl_sf_exp(-Xbetai)); } } - // The gradient of f, df = (df/dx, df/dy). -void wgsl_mixed_optim_df (const gsl_vector *beta, void *params, - gsl_vector *out) { +void wgsl_mixed_optim_df(const gsl_vector *beta, void *params, + gsl_vector *out) { fix_parm_mixed_T *p = (fix_parm_mixed_T *)params; int n = p->y->size; int K = p->X->size2; @@ -91,50 +87,49 @@ void wgsl_mixed_optim_df (const gsl_vector *beta, void *params, int npar = beta->size; // Intitialize gradient out necessary? - for(int i = 0; i < npar; ++i) - out->data[i]= 0; + for (int i = 0; i < npar; ++i) + out->data[i] = 0; // Changed loop start at 1 instead of 0 to avoid regularization of beta 0. - for(int i = 1; i < npar; ++i) - out->data[i]= p->lambdaL2*beta->data[i]; - for(int i = 1; i < npar; ++i) - out->data[i]+= p->lambdaL1*((beta->data[i]>0)-(beta->data[i]<0)); - - for(int i = 0; i < n; ++i) { - double pn=0; - double Xbetai=beta->data[0]; - int iParm=1; - for(int k = 0; k < K; ++k) { - if(gsl_matrix_int_get(p->X,i,k)>0) - Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm]; - iParm+=p->nlev->data[k]-1; + for (int i = 1; i < npar; ++i) + out->data[i] = p->lambdaL2 * beta->data[i]; + for (int i = 1; i < npar; ++i) + out->data[i] += p->lambdaL1 * ((beta->data[i] > 0) - (beta->data[i] < 0)); + + for (int i = 0; i < n; ++i) { + double pn = 0; + double Xbetai = beta->data[0]; + int iParm = 1; + for (int k = 0; k < K; ++k) { + if (gsl_matrix_int_get(p->X, i, k) > 0) + Xbetai += beta->data[gsl_matrix_int_get(p->X, i, k) - 1 + iParm]; + iParm += p->nlev->data[k] - 1; } // Adding the continuous. - for(int k = 0; k < Kc; ++k) - Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm++]; + for (int k = 0; k < Kc; ++k) + Xbetai += gsl_matrix_get(p->Xc, i, k) * beta->data[iParm++]; - pn= -( p->y->data[i] - 1/(1 + gsl_sf_exp(-Xbetai)) ); + pn = -(p->y->data[i] - 1 / (1 + gsl_sf_exp(-Xbetai))); - out->data[0]+= pn; - iParm=1; - for(int k = 0; k < K; ++k) { - if(gsl_matrix_int_get(p->X,i,k)>0) - out->data[gsl_matrix_int_get(p->X,i,k)-1+iParm]+=pn; - iParm+=p->nlev->data[k]-1; + out->data[0] += pn; + iParm = 1; + for (int k = 0; k < K; ++k) { + if (gsl_matrix_int_get(p->X, i, k) > 0) + out->data[gsl_matrix_int_get(p->X, i, k) - 1 + iParm] += pn; + iParm += p->nlev->data[k] - 1; } // Adding the continuous. - for(int k = 0; k < Kc; ++k) { - out->data[iParm++] += gsl_matrix_get(p->Xc,i,k)*pn; + for (int k = 0; k < Kc; ++k) { + out->data[iParm++] += gsl_matrix_get(p->Xc, i, k) * pn; } } - } // The Hessian of f. -void wgsl_mixed_optim_hessian (const gsl_vector *beta, void *params, - gsl_matrix *out) { +void wgsl_mixed_optim_hessian(const gsl_vector *beta, void *params, + gsl_matrix *out) { fix_parm_mixed_T *p = (fix_parm_mixed_T *)params; int n = p->y->size; int K = p->X->size2; @@ -146,120 +141,121 @@ void wgsl_mixed_optim_hessian (const gsl_vector *beta, void *params, gsl_matrix_set_zero(out); /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0*/ - for(int i = 1; i < npar; ++i) - gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double-check this. + for (int i = 1; i < npar; ++i) + gsl_matrix_set(out, i, i, (p->lambdaL2)); // Double-check this. // L1 penalty not working yet, as not differentiable, I may need to // do coordinate descent (as in glm_net) - for(int i = 0; i < n; ++i) { - double pn=0; - double aux=0; - double Xbetai=beta->data[0]; - int iParm1=1; - for(int k = 0; k < K; ++k) { - if(gsl_matrix_int_get(p->X,i,k)>0) - Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm1]; - iParm1+=p->nlev->data[k]-1; //-1? + for (int i = 0; i < n; ++i) { + double pn = 0; + double aux = 0; + double Xbetai = beta->data[0]; + int iParm1 = 1; + for (int k = 0; k < K; ++k) { + if (gsl_matrix_int_get(p->X, i, k) > 0) + Xbetai += beta->data[gsl_matrix_int_get(p->X, i, k) - 1 + iParm1]; + iParm1 += p->nlev->data[k] - 1; //-1? } // Adding the continuous. - for(int k = 0; k < Kc; ++k) - Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm1++]; + for (int k = 0; k < Kc; ++k) + Xbetai += gsl_matrix_get(p->Xc, i, k) * beta->data[iParm1++]; - pn= 1/(1 + gsl_sf_exp(-Xbetai)); + pn = 1 / (1 + gsl_sf_exp(-Xbetai)); // Add a protection for pn very close to 0 or 1? - aux=pn*(1-pn); + aux = pn * (1 - pn); // Calculate sub-gradient vector gn. gsl_vector_set_zero(gn); - gn->data[0]= 1; - iParm1=1; - for(int k = 0; k < K; ++k) { - if(gsl_matrix_int_get(p->X,i,k)>0) - gn->data[gsl_matrix_int_get(p->X,i,k)-1+iParm1]=1; - iParm1+=p->nlev->data[k]-1; + gn->data[0] = 1; + iParm1 = 1; + for (int k = 0; k < K; ++k) { + if (gsl_matrix_int_get(p->X, i, k) > 0) + gn->data[gsl_matrix_int_get(p->X, i, k) - 1 + iParm1] = 1; + iParm1 += p->nlev->data[k] - 1; } // Adding the continuous. - for(int k = 0; k < Kc; ++k) { - gn->data[iParm1++] = gsl_matrix_get(p->Xc,i,k); + for (int k = 0; k < Kc; ++k) { + gn->data[iParm1++] = gsl_matrix_get(p->Xc, i, k); } - for(int k1=0;k1<npar; ++k1) - if(gn->data[k1]!=0) - for(int k2=0;k2<npar; ++k2) - if(gn->data[k2]!=0) - *gsl_matrix_ptr(out,k1,k2) += (aux * gn->data[k1] * gn->data[k2]); + for (int k1 = 0; k1 < npar; ++k1) + if (gn->data[k1] != 0) + for (int k2 = 0; k2 < npar; ++k2) + if (gn->data[k2] != 0) + *gsl_matrix_ptr(out, k1, k2) += (aux * gn->data[k1] * gn->data[k2]); } gsl_vector_free(gn); } double wgsl_mixed_optim_f(gsl_vector *v, void *params) { - double mLogLik=0; + double mLogLik = 0; fix_parm_mixed_T *p = (fix_parm_mixed_T *)params; - mLogLik = fLogit_mixed(v,p->X,p->nlev,p->Xc,p->y,p->lambdaL1,p->lambdaL2); + mLogLik = + fLogit_mixed(v, p->X, p->nlev, p->Xc, p->y, p->lambdaL1, p->lambdaL2); return mLogLik; } // Compute both f and df together. -void -wgsl_mixed_optim_fdf (gsl_vector *x, void *params, double *f, gsl_vector *df) { +void wgsl_mixed_optim_fdf(gsl_vector *x, void *params, double *f, + gsl_vector *df) { *f = wgsl_mixed_optim_f(x, params); wgsl_mixed_optim_df(x, params, df); } // Xc is the matrix of continuous covariates, Nobs x Kc (NULL if not used). int logistic_mixed_fit(gsl_vector *beta, gsl_matrix_int *X, - gsl_vector_int *nlev, gsl_matrix *Xc, - gsl_vector *y, double lambdaL1, double lambdaL2) { - double mLogLik=0; + gsl_vector_int *nlev, gsl_matrix *Xc, gsl_vector *y, + double lambdaL1, double lambdaL2) { + double mLogLik = 0; fix_parm_mixed_T p; int npar = beta->size; - int iter=0; - double maxchange=0; + int iter = 0; + double maxchange = 0; // Intializing fix parameters. - p.X=X; - p.Xc=Xc; - p.nlev=nlev; - p.y=y; - p.lambdaL1=lambdaL1; - p.lambdaL2=lambdaL2; + p.X = X; + p.Xc = Xc; + p.nlev = nlev; + p.y = y; + p.lambdaL1 = lambdaL1; + p.lambdaL2 = lambdaL2; // Initial fit. - mLogLik = wgsl_mixed_optim_f(beta,&p); + mLogLik = wgsl_mixed_optim_f(beta, &p); - gsl_matrix *myH = gsl_matrix_alloc(npar,npar); // Hessian matrix. - gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. + gsl_matrix *myH = gsl_matrix_alloc(npar, npar); // Hessian matrix. + gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. - gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. - gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. + gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. + gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. - for(iter=0;iter<100;iter++){ - wgsl_mixed_optim_hessian(beta,&p,myH); // Calculate Hessian. - wgsl_mixed_optim_df(beta,&p,myG); // Calculate Gradient. - gsl_linalg_QR_decomp(myH,tau); // Calculate next beta. - gsl_linalg_QR_solve(myH,tau,myG,stBeta); - gsl_vector_sub(beta,stBeta); + for (iter = 0; iter < 100; iter++) { + wgsl_mixed_optim_hessian(beta, &p, myH); // Calculate Hessian. + wgsl_mixed_optim_df(beta, &p, myG); // Calculate Gradient. + gsl_linalg_QR_decomp(myH, tau); // Calculate next beta. + gsl_linalg_QR_solve(myH, tau, myG, stBeta); + gsl_vector_sub(beta, stBeta); // Monitor convergence. - maxchange=0; - for(int i=0;i<npar; i++) - if(maxchange<fabs(stBeta->data[i])) - maxchange=fabs(stBeta->data[i]); + maxchange = 0; + for (int i = 0; i < npar; i++) + if (maxchange < fabs(stBeta->data[i])) + maxchange = fabs(stBeta->data[i]); - if(maxchange<1E-4) + if (maxchange < 1E-4) break; } // Final fit. - mLogLik = wgsl_mixed_optim_f(beta,&p); + mLogLik = wgsl_mixed_optim_f(beta, &p); - gsl_vector_free (tau); - gsl_vector_free (stBeta); - gsl_vector_free (myG); - gsl_matrix_free (myH); + gsl_vector_free(tau); + gsl_vector_free(stBeta); + gsl_vector_free(myG); + gsl_matrix_free(myH); return 0; } @@ -278,8 +274,8 @@ typedef struct { double lambdaL2; } fix_parm_cat_T; -double fLogit_cat (gsl_vector *beta, gsl_matrix_int *X, gsl_vector_int *nlev, - gsl_vector *y, double lambdaL1, double lambdaL2) { +double fLogit_cat(gsl_vector *beta, gsl_matrix_int *X, gsl_vector_int *nlev, + gsl_vector *y, double lambdaL1, double lambdaL2) { int n = y->size; int npar = beta->size; double total = 0; @@ -288,91 +284,90 @@ double fLogit_cat (gsl_vector *beta, gsl_matrix_int *X, gsl_vector_int *nlev, // omp_set_num_threads(ompthr); /\* Changed loop start at 1 instead // of 0 to avoid regularization of beta 0*\/ /\*#pragma omp parallel // for reduction (+:total)*\/ - for(int i = 1; i < npar; ++i) - total += beta->data[i]*beta->data[i]; - total = (-total*lambdaL2/2); + for (int i = 1; i < npar; ++i) + total += beta->data[i] * beta->data[i]; + total = (-total * lambdaL2 / 2); // /\*#pragma omp parallel for reduction (+:aux)*\/ - for(int i = 1; i < npar; ++i) - aux += (beta->data[i]>0 ? beta->data[i] : -beta->data[i]); - total = total-aux*lambdaL1; + for (int i = 1; i < npar; ++i) + aux += (beta->data[i] > 0 ? beta->data[i] : -beta->data[i]); + total = total - aux * lambdaL1; // #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y) // #reduction (+:total) - for(int i = 0; i < n; ++i) { - double Xbetai=beta->data[0]; - int iParm=1; - for(int k = 0; k < X->size2; ++k) { - if(gsl_matrix_int_get(X,i,k)>0) - Xbetai+=beta->data[gsl_matrix_int_get(X,i,k)-1+iParm]; - iParm+=nlev->data[k]-1; + for (int i = 0; i < n; ++i) { + double Xbetai = beta->data[0]; + int iParm = 1; + for (int k = 0; k < X->size2; ++k) { + if (gsl_matrix_int_get(X, i, k) > 0) + Xbetai += beta->data[gsl_matrix_int_get(X, i, k) - 1 + iParm]; + iParm += nlev->data[k] - 1; } - total += y->data[i]*Xbetai-gsl_sf_log_1plusx(gsl_sf_exp(Xbetai)); + total += y->data[i] * Xbetai - gsl_sf_log_1plusx(gsl_sf_exp(Xbetai)); } return -total; } -void logistic_cat_pred (gsl_vector *beta, // Vector of parameters - // length = 1 + Sum_k(C_k-1). - gsl_matrix_int *X, // Matrix Nobs x K - gsl_vector_int *nlev, // Vector with #categories - gsl_vector *yhat){ // Vector of prob. predicted by - // the logistic. - for(int i = 0; i < X->size1; ++i) { - double Xbetai=beta->data[0]; - int iParm=1; - for(int k = 0; k < X->size2; ++k) { - if(gsl_matrix_int_get(X,i,k)>0) - Xbetai+=beta->data[gsl_matrix_int_get(X,i,k)-1+iParm]; - iParm+=nlev->data[k]-1; +void logistic_cat_pred(gsl_vector *beta, // Vector of parameters + // length = 1 + Sum_k(C_k-1). + gsl_matrix_int *X, // Matrix Nobs x K + gsl_vector_int *nlev, // Vector with #categories + gsl_vector *yhat) { // Vector of prob. predicted by + // the logistic. + for (int i = 0; i < X->size1; ++i) { + double Xbetai = beta->data[0]; + int iParm = 1; + for (int k = 0; k < X->size2; ++k) { + if (gsl_matrix_int_get(X, i, k) > 0) + Xbetai += beta->data[gsl_matrix_int_get(X, i, k) - 1 + iParm]; + iParm += nlev->data[k] - 1; } - yhat->data[i]=1/(1 + gsl_sf_exp(-Xbetai)); + yhat->data[i] = 1 / (1 + gsl_sf_exp(-Xbetai)); } } // The gradient of f, df = (df/dx, df/dy). -void wgsl_cat_optim_df (const gsl_vector *beta, void *params, - gsl_vector *out) { +void wgsl_cat_optim_df(const gsl_vector *beta, void *params, gsl_vector *out) { fix_parm_cat_T *p = (fix_parm_cat_T *)params; int n = p->y->size; int K = p->X->size2; int npar = beta->size; // Intitialize gradient out necessary? - for(int i = 0; i < npar; ++i) - out->data[i]= 0; + for (int i = 0; i < npar; ++i) + out->data[i] = 0; // Changed loop start at 1 instead of 0 to avoid regularization of beta 0. - for(int i = 1; i < npar; ++i) - out->data[i]= p->lambdaL2*beta->data[i]; - for(int i = 1; i < npar; ++i) - out->data[i]+= p->lambdaL1*((beta->data[i]>0)-(beta->data[i]<0)); - - for(int i = 0; i < n; ++i) { - double pn=0; - double Xbetai=beta->data[0]; - int iParm=1; - for(int k = 0; k < K; ++k) { - if(gsl_matrix_int_get(p->X,i,k)>0) - Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm]; - iParm+=p->nlev->data[k]-1; + for (int i = 1; i < npar; ++i) + out->data[i] = p->lambdaL2 * beta->data[i]; + for (int i = 1; i < npar; ++i) + out->data[i] += p->lambdaL1 * ((beta->data[i] > 0) - (beta->data[i] < 0)); + + for (int i = 0; i < n; ++i) { + double pn = 0; + double Xbetai = beta->data[0]; + int iParm = 1; + for (int k = 0; k < K; ++k) { + if (gsl_matrix_int_get(p->X, i, k) > 0) + Xbetai += beta->data[gsl_matrix_int_get(p->X, i, k) - 1 + iParm]; + iParm += p->nlev->data[k] - 1; } - pn= -( p->y->data[i] - 1/(1 + gsl_sf_exp(-Xbetai)) ); + pn = -(p->y->data[i] - 1 / (1 + gsl_sf_exp(-Xbetai))); - out->data[0]+= pn; - iParm=1; - for(int k = 0; k < K; ++k) { - if(gsl_matrix_int_get(p->X,i,k)>0) - out->data[gsl_matrix_int_get(p->X,i,k)-1+iParm]+=pn; - iParm+=p->nlev->data[k]-1; + out->data[0] += pn; + iParm = 1; + for (int k = 0; k < K; ++k) { + if (gsl_matrix_int_get(p->X, i, k) > 0) + out->data[gsl_matrix_int_get(p->X, i, k) - 1 + iParm] += pn; + iParm += p->nlev->data[k] - 1; } } } // The Hessian of f. -void wgsl_cat_optim_hessian (const gsl_vector *beta, void *params, - gsl_matrix *out) { +void wgsl_cat_optim_hessian(const gsl_vector *beta, void *params, + gsl_matrix *out) { fix_parm_cat_T *p = (fix_parm_cat_T *)params; int n = p->y->size; int K = p->X->size2; @@ -382,123 +377,119 @@ void wgsl_cat_optim_hessian (const gsl_vector *beta, void *params, gsl_matrix_set_zero(out); // Changed loop start at 1 instead of 0 to avoid regularization of beta. - for(int i = 1; i < npar; ++i) - gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double-check this. + for (int i = 1; i < npar; ++i) + gsl_matrix_set(out, i, i, (p->lambdaL2)); // Double-check this. // L1 penalty not working yet, as not differentiable, I may need to // do coordinate descent (as in glm_net). - for(int i = 0; i < n; ++i) { - double pn=0; - double aux=0; - double Xbetai=beta->data[0]; - int iParm2=1; - int iParm1=1; - for(int k = 0; k < K; ++k) { - if(gsl_matrix_int_get(p->X,i,k)>0) - Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm1]; - iParm1+=p->nlev->data[k]-1; //-1? + for (int i = 0; i < n; ++i) { + double pn = 0; + double aux = 0; + double Xbetai = beta->data[0]; + int iParm2 = 1; + int iParm1 = 1; + for (int k = 0; k < K; ++k) { + if (gsl_matrix_int_get(p->X, i, k) > 0) + Xbetai += beta->data[gsl_matrix_int_get(p->X, i, k) - 1 + iParm1]; + iParm1 += p->nlev->data[k] - 1; //-1? } - pn= 1/(1 + gsl_sf_exp(-Xbetai)); + pn = 1 / (1 + gsl_sf_exp(-Xbetai)); // Add a protection for pn very close to 0 or 1? - aux=pn*(1-pn); - *gsl_matrix_ptr(out,0,0)+=aux; - iParm2=1; - for(int k2 = 0; k2 < K; ++k2) { - if(gsl_matrix_int_get(p->X,i,k2)>0) - *gsl_matrix_ptr(out,0,gsl_matrix_int_get(p->X,i,k2)-1+iParm2)+=aux; - iParm2+=p->nlev->data[k2]-1; //-1? + aux = pn * (1 - pn); + *gsl_matrix_ptr(out, 0, 0) += aux; + iParm2 = 1; + for (int k2 = 0; k2 < K; ++k2) { + if (gsl_matrix_int_get(p->X, i, k2) > 0) + *gsl_matrix_ptr(out, 0, gsl_matrix_int_get(p->X, i, k2) - 1 + iParm2) += + aux; + iParm2 += p->nlev->data[k2] - 1; //-1? } - iParm1=1; - for(int k1 = 0; k1 < K; ++k1) { - if(gsl_matrix_int_get(p->X,i,k1)>0) - *gsl_matrix_ptr(out,gsl_matrix_int_get(p->X,i,k1)-1+iParm1,0)+=aux; - iParm2=1; - for(int k2 = 0; k2 < K; ++k2) { - if((gsl_matrix_int_get(p->X,i,k1)>0) && - (gsl_matrix_int_get(p->X,i,k2)>0)) - *gsl_matrix_ptr(out - ,gsl_matrix_int_get(p->X,i,k1)-1+iParm1 - ,gsl_matrix_int_get(p->X,i,k2)-1+iParm2 - )+=aux; - iParm2+=p->nlev->data[k2]-1; //-1? + iParm1 = 1; + for (int k1 = 0; k1 < K; ++k1) { + if (gsl_matrix_int_get(p->X, i, k1) > 0) + *gsl_matrix_ptr(out, gsl_matrix_int_get(p->X, i, k1) - 1 + iParm1, 0) += + aux; + iParm2 = 1; + for (int k2 = 0; k2 < K; ++k2) { + if ((gsl_matrix_int_get(p->X, i, k1) > 0) && + (gsl_matrix_int_get(p->X, i, k2) > 0)) + *gsl_matrix_ptr(out, gsl_matrix_int_get(p->X, i, k1) - 1 + iParm1, + gsl_matrix_int_get(p->X, i, k2) - 1 + iParm2) += aux; + iParm2 += p->nlev->data[k2] - 1; //-1? } - iParm1+=p->nlev->data[k1]-1; //-1? + iParm1 += p->nlev->data[k1] - 1; //-1? } } } double wgsl_cat_optim_f(gsl_vector *v, void *params) { - double mLogLik=0; + double mLogLik = 0; fix_parm_cat_T *p = (fix_parm_cat_T *)params; - mLogLik = fLogit_cat(v,p->X,p->nlev,p->y,p->lambdaL1,p->lambdaL2); + mLogLik = fLogit_cat(v, p->X, p->nlev, p->y, p->lambdaL1, p->lambdaL2); return mLogLik; } // Compute both f and df together. -void wgsl_cat_optim_fdf (gsl_vector *x, void *params, double *f, - gsl_vector *df) { +void wgsl_cat_optim_fdf(gsl_vector *x, void *params, double *f, + gsl_vector *df) { *f = wgsl_cat_optim_f(x, params); wgsl_cat_optim_df(x, params, df); } -int logistic_cat_fit(gsl_vector *beta, - gsl_matrix_int *X, - gsl_vector_int *nlev, - gsl_vector *y, - double lambdaL1, - double lambdaL2) { - double mLogLik=0; +int logistic_cat_fit(gsl_vector *beta, gsl_matrix_int *X, gsl_vector_int *nlev, + gsl_vector *y, double lambdaL1, double lambdaL2) { + double mLogLik = 0; fix_parm_cat_T p; int npar = beta->size; - int iter=0; - double maxchange=0; + int iter = 0; + double maxchange = 0; // Intializing fix parameters. - p.X=X; - p.nlev=nlev; - p.y=y; - p.lambdaL1=lambdaL1; - p.lambdaL2=lambdaL2; + p.X = X; + p.nlev = nlev; + p.y = y; + p.lambdaL1 = lambdaL1; + p.lambdaL2 = lambdaL2; // Initial fit. - mLogLik = wgsl_cat_optim_f(beta,&p); + mLogLik = wgsl_cat_optim_f(beta, &p); - gsl_matrix *myH = gsl_matrix_alloc(npar,npar); // Hessian matrix. - gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. + gsl_matrix *myH = gsl_matrix_alloc(npar, npar); // Hessian matrix. + gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. - gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. - gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. + gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. + gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. - for(iter=0;iter<100;iter++){ - wgsl_cat_optim_hessian(beta,&p,myH); // Calculate Hessian. - wgsl_cat_optim_df(beta,&p,myG); // Calculate Gradient. - gsl_linalg_QR_decomp(myH,tau); // Calculate next beta. - gsl_linalg_QR_solve(myH,tau,myG,stBeta); - gsl_vector_sub(beta,stBeta); + for (iter = 0; iter < 100; iter++) { + wgsl_cat_optim_hessian(beta, &p, myH); // Calculate Hessian. + wgsl_cat_optim_df(beta, &p, myG); // Calculate Gradient. + gsl_linalg_QR_decomp(myH, tau); // Calculate next beta. + gsl_linalg_QR_solve(myH, tau, myG, stBeta); + gsl_vector_sub(beta, stBeta); // Monitor convergence. - maxchange=0; - for(int i=0;i<npar; i++) - if(maxchange<fabs(stBeta->data[i])) - maxchange=fabs(stBeta->data[i]); + maxchange = 0; + for (int i = 0; i < npar; i++) + if (maxchange < fabs(stBeta->data[i])) + maxchange = fabs(stBeta->data[i]); #ifdef _RPR_DEBUG_ - mLogLik = wgsl_cat_optim_f(beta,&p); + mLogLik = wgsl_cat_optim_f(beta, &p); #endif - if(maxchange<1E-4) + if (maxchange < 1E-4) break; } // Final fit. - mLogLik = wgsl_cat_optim_f(beta,&p); + mLogLik = wgsl_cat_optim_f(beta, &p); - gsl_vector_free (tau); - gsl_vector_free (stBeta); - gsl_vector_free (myG); - gsl_matrix_free (myH); + gsl_vector_free(tau); + gsl_vector_free(stBeta); + gsl_vector_free(myG); + gsl_matrix_free(myH); return 0; } @@ -509,15 +500,15 @@ int logistic_cat_fit(gsl_vector *beta, // I need to bundle all the data that goes to the function to optimze // together. -typedef struct{ - gsl_matrix *Xc; // continuous covariates; Matrix Nobs x Kc +typedef struct { + gsl_matrix *Xc; // continuous covariates; Matrix Nobs x Kc gsl_vector *y; double lambdaL1; double lambdaL2; -}fix_parm_cont_T; +} fix_parm_cont_T; double fLogit_cont(gsl_vector *beta, gsl_matrix *Xc, gsl_vector *y, - double lambdaL1, double lambdaL2) { + double lambdaL1, double lambdaL2) { int n = y->size; int npar = beta->size; double total = 0; @@ -526,82 +517,81 @@ double fLogit_cont(gsl_vector *beta, gsl_matrix *Xc, gsl_vector *y, // omp_set_num_threads(ompthr); /\* Changed loop start at 1 instead // of 0 to avoid regularization of beta_0*\/ /\*#pragma omp parallel // for reduction (+:total)*\/ - for(int i = 1; i < npar; ++i) - total += beta->data[i]*beta->data[i]; - total = (-total*lambdaL2/2); + for (int i = 1; i < npar; ++i) + total += beta->data[i] * beta->data[i]; + total = (-total * lambdaL2 / 2); // /\*#pragma omp parallel for reduction (+:aux)*\/ - for(int i = 1; i < npar; ++i) - aux += (beta->data[i]>0 ? beta->data[i] : -beta->data[i]); - total = total-aux*lambdaL1; + for (int i = 1; i < npar; ++i) + aux += (beta->data[i] > 0 ? beta->data[i] : -beta->data[i]); + total = total - aux * lambdaL1; // #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y) // #reduction (+:total) - for(int i = 0; i < n; ++i) { - double Xbetai=beta->data[0]; - int iParm=1; - for(int k = 0; k < (Xc->size2); ++k) - Xbetai+= gsl_matrix_get(Xc,i,k)*beta->data[iParm++]; - total += y->data[i]*Xbetai-gsl_sf_log_1plusx(gsl_sf_exp(Xbetai)); + for (int i = 0; i < n; ++i) { + double Xbetai = beta->data[0]; + int iParm = 1; + for (int k = 0; k < (Xc->size2); ++k) + Xbetai += gsl_matrix_get(Xc, i, k) * beta->data[iParm++]; + total += y->data[i] * Xbetai - gsl_sf_log_1plusx(gsl_sf_exp(Xbetai)); } return -total; } -void logistic_cont_pred(gsl_vector *beta, // Vector of parameters - // length = 1 + Sum_k(C_k-1). - gsl_matrix *Xc, // Continuous covariates matrix, - // Nobs x Kc (NULL if not used). - gsl_vector *yhat) { // Vector of prob. predicted by - // the logistic. - for(int i = 0; i < Xc->size1; ++i) { - double Xbetai=beta->data[0]; - int iParm=1; - for(int k = 0; k < (Xc->size2); ++k) - Xbetai+= gsl_matrix_get(Xc,i,k)*beta->data[iParm++]; - yhat->data[i]=1/(1 + gsl_sf_exp(-Xbetai)); +void logistic_cont_pred(gsl_vector *beta, // Vector of parameters + // length = 1 + Sum_k(C_k-1). + gsl_matrix *Xc, // Continuous covariates matrix, + // Nobs x Kc (NULL if not used). + gsl_vector *yhat) { // Vector of prob. predicted by + // the logistic. + for (int i = 0; i < Xc->size1; ++i) { + double Xbetai = beta->data[0]; + int iParm = 1; + for (int k = 0; k < (Xc->size2); ++k) + Xbetai += gsl_matrix_get(Xc, i, k) * beta->data[iParm++]; + yhat->data[i] = 1 / (1 + gsl_sf_exp(-Xbetai)); } } // The gradient of f, df = (df/dx, df/dy). -void wgsl_cont_optim_df (const gsl_vector *beta, void *params, - gsl_vector *out) { +void wgsl_cont_optim_df(const gsl_vector *beta, void *params, gsl_vector *out) { fix_parm_cont_T *p = (fix_parm_cont_T *)params; int n = p->y->size; int Kc = p->Xc->size2; int npar = beta->size; // Intitialize gradient out necessary? - for(int i = 0; i < npar; ++i) - out->data[i]= 0; + for (int i = 0; i < npar; ++i) + out->data[i] = 0; // Changed loop start at 1 instead of 0 to avoid regularization of beta 0. - for(int i = 1; i < npar; ++i) - out->data[i]= p->lambdaL2*beta->data[i]; - for(int i = 1; i < npar; ++i) - out->data[i]+= p->lambdaL1*((beta->data[i]>0)-(beta->data[i]<0)); + for (int i = 1; i < npar; ++i) + out->data[i] = p->lambdaL2 * beta->data[i]; + for (int i = 1; i < npar; ++i) + out->data[i] += p->lambdaL1 * ((beta->data[i] > 0) - (beta->data[i] < 0)); - for(int i = 0; i < n; ++i) { - double pn=0; - double Xbetai=beta->data[0]; - int iParm=1; - for(int k = 0; k < Kc; ++k) - Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm++]; + for (int i = 0; i < n; ++i) { + double pn = 0; + double Xbetai = beta->data[0]; + int iParm = 1; + for (int k = 0; k < Kc; ++k) + Xbetai += gsl_matrix_get(p->Xc, i, k) * beta->data[iParm++]; - pn= -( p->y->data[i] - 1/(1 + gsl_sf_exp(-Xbetai)) ); + pn = -(p->y->data[i] - 1 / (1 + gsl_sf_exp(-Xbetai))); - out->data[0]+= pn; - iParm=1; + out->data[0] += pn; + iParm = 1; // Adding the continuous. - for(int k = 0; k < Kc; ++k) { - out->data[iParm++] += gsl_matrix_get(p->Xc,i,k)*pn; + for (int k = 0; k < Kc; ++k) { + out->data[iParm++] += gsl_matrix_get(p->Xc, i, k) * pn; } } } // The Hessian of f. -void wgsl_cont_optim_hessian (const gsl_vector *beta, void *params, - gsl_matrix *out) { +void wgsl_cont_optim_hessian(const gsl_vector *beta, void *params, + gsl_matrix *out) { fix_parm_cont_T *p = (fix_parm_cont_T *)params; int n = p->y->size; int Kc = p->Xc->size2; @@ -614,111 +604,109 @@ void wgsl_cont_optim_hessian (const gsl_vector *beta, void *params, // Changed loop start at 1 instead of 0 to avoid regularization of // beta 0. - for(int i = 1; i < npar; ++i) - gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double-check this. + for (int i = 1; i < npar; ++i) + gsl_matrix_set(out, i, i, (p->lambdaL2)); // Double-check this. // L1 penalty not working yet, as not differentiable, I may need to // do coordinate descent (as in glm_net). - for(int i = 0; i < n; ++i) { - double pn=0; - double aux=0; - double Xbetai=beta->data[0]; - int iParm1=1; - for(int k = 0; k < Kc; ++k) - Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm1++]; + for (int i = 0; i < n; ++i) { + double pn = 0; + double aux = 0; + double Xbetai = beta->data[0]; + int iParm1 = 1; + for (int k = 0; k < Kc; ++k) + Xbetai += gsl_matrix_get(p->Xc, i, k) * beta->data[iParm1++]; - pn= 1/(1 + gsl_sf_exp(-Xbetai)); + pn = 1 / (1 + gsl_sf_exp(-Xbetai)); // Add a protection for pn very close to 0 or 1? - aux=pn*(1-pn); + aux = pn * (1 - pn); // Calculate sub-gradient vector gn. gsl_vector_set_zero(gn); - gn->data[0]= 1; - iParm1=1; - for(int k = 0; k < Kc; ++k) { - gn->data[iParm1++] = gsl_matrix_get(p->Xc,i,k); + gn->data[0] = 1; + iParm1 = 1; + for (int k = 0; k < Kc; ++k) { + gn->data[iParm1++] = gsl_matrix_get(p->Xc, i, k); } - for(int k1=0;k1<npar; ++k1) - if(gn->data[k1]!=0) - for(int k2=0;k2<npar; ++k2) - if(gn->data[k2]!=0) - *gsl_matrix_ptr(out,k1,k2) += (aux * gn->data[k1] * gn->data[k2]); + for (int k1 = 0; k1 < npar; ++k1) + if (gn->data[k1] != 0) + for (int k2 = 0; k2 < npar; ++k2) + if (gn->data[k2] != 0) + *gsl_matrix_ptr(out, k1, k2) += (aux * gn->data[k1] * gn->data[k2]); } gsl_vector_free(gn); } double wgsl_cont_optim_f(gsl_vector *v, void *params) { - double mLogLik=0; + double mLogLik = 0; fix_parm_cont_T *p = (fix_parm_cont_T *)params; - mLogLik = fLogit_cont(v,p->Xc,p->y,p->lambdaL1,p->lambdaL2); + mLogLik = fLogit_cont(v, p->Xc, p->y, p->lambdaL1, p->lambdaL2); return mLogLik; } // Compute both f and df together. -void wgsl_cont_optim_fdf (gsl_vector *x, void *params, - double *f, gsl_vector *df) { +void wgsl_cont_optim_fdf(gsl_vector *x, void *params, double *f, + gsl_vector *df) { *f = wgsl_cont_optim_f(x, params); wgsl_cont_optim_df(x, params, df); } -int logistic_cont_fit (gsl_vector *beta, - gsl_matrix *Xc, // Continuous covariates matrix, - // Nobs x Kc (NULL if not used). - gsl_vector *y, - double lambdaL1, - double lambdaL2) { +int logistic_cont_fit(gsl_vector *beta, + gsl_matrix *Xc, // Continuous covariates matrix, + // Nobs x Kc (NULL if not used). + gsl_vector *y, double lambdaL1, double lambdaL2) { - double mLogLik=0; + double mLogLik = 0; fix_parm_cont_T p; int npar = beta->size; - int iter=0; - double maxchange=0; + int iter = 0; + double maxchange = 0; // Initializing fix parameters. - p.Xc=Xc; - p.y=y; - p.lambdaL1=lambdaL1; - p.lambdaL2=lambdaL2; + p.Xc = Xc; + p.y = y; + p.lambdaL1 = lambdaL1; + p.lambdaL2 = lambdaL2; // Initial fit. - mLogLik = wgsl_cont_optim_f(beta,&p); + mLogLik = wgsl_cont_optim_f(beta, &p); - gsl_matrix *myH = gsl_matrix_alloc(npar,npar); // Hessian matrix. - gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. + gsl_matrix *myH = gsl_matrix_alloc(npar, npar); // Hessian matrix. + gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. - gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. - gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. + gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. + gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. - for(iter=0;iter<100;iter++){ - wgsl_cont_optim_hessian(beta,&p,myH); // Calculate Hessian. - wgsl_cont_optim_df(beta,&p,myG); // Calculate Gradient. - gsl_linalg_QR_decomp(myH,tau); // Calculate next beta. - gsl_linalg_QR_solve(myH,tau,myG,stBeta); - gsl_vector_sub(beta,stBeta); + for (iter = 0; iter < 100; iter++) { + wgsl_cont_optim_hessian(beta, &p, myH); // Calculate Hessian. + wgsl_cont_optim_df(beta, &p, myG); // Calculate Gradient. + gsl_linalg_QR_decomp(myH, tau); // Calculate next beta. + gsl_linalg_QR_solve(myH, tau, myG, stBeta); + gsl_vector_sub(beta, stBeta); // Monitor convergence. - maxchange=0; - for(int i=0;i<npar; i++) - if(maxchange<fabs(stBeta->data[i])) - maxchange=fabs(stBeta->data[i]); + maxchange = 0; + for (int i = 0; i < npar; i++) + if (maxchange < fabs(stBeta->data[i])) + maxchange = fabs(stBeta->data[i]); #ifdef _RPR_DEBUG_ - mLogLik = wgsl_cont_optim_f(beta,&p); + mLogLik = wgsl_cont_optim_f(beta, &p); #endif - if(maxchange<1E-4) + if (maxchange < 1E-4) break; } // Final fit. - mLogLik = wgsl_cont_optim_f(beta,&p); + mLogLik = wgsl_cont_optim_f(beta, &p); - gsl_vector_free (tau); - gsl_vector_free (stBeta); - gsl_vector_free (myG); - gsl_matrix_free (myH); + gsl_vector_free(tau); + gsl_vector_free(stBeta); + gsl_vector_free(myG); + gsl_matrix_free(myH); return 0; } |