diff options
Diffstat (limited to 'src/logistic.cpp')
| -rw-r--r-- | src/logistic.cpp | 783 |
1 files changed, 783 insertions, 0 deletions
diff --git a/src/logistic.cpp b/src/logistic.cpp new file mode 100644 index 0000000..1b47946 --- /dev/null +++ b/src/logistic.cpp @@ -0,0 +1,783 @@ +#include <stdio.h>
+#include <math.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_linalg.h>
+#include "logistic.h"
+
+// I need to bundle all the data that goes to the function to optimze together.
+typedef struct{
+ gsl_matrix_int *X;
+ gsl_vector_int *nlev;
+ gsl_vector *y;
+ gsl_matrix *Xc; // continuous covariates Matrix Nobs x Kc (NULL if not used)
+ double lambdaL1;
+ 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)
+{
+ int n = y->size;
+ // int k = X->size2;
+ int npar = beta->size;
+ double total = 0;
+ double aux = 0;
+ /* 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);
+ /* /\*#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;
+ /* #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 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 Nobs 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));
+ }
+}
+
+
+/* The gradient of f, df = (df/dx, df/dy). */
+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;
+ int Kc = p->Xc->size2;
+ int npar = beta->size;
+ // Intitialize gradient out necessary?
+ 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;
+ }
+ // Adding the continuous
+ 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)) );
+
+ 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;
+ }
+ }
+
+}
+
+
+/* The Hessian of f */
+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;
+ int Kc = p->Xc->size2;
+ int npar = beta->size;
+ gsl_vector *gn = gsl_vector_alloc(npar); /* gn */
+ // Intitialize Hessian out necessary ???
+ 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
+ // 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?
+ }
+ // Adding the continuous
+ 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));
+ // Add a protection for pn very close to 0 or 1?
+ 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;
+ }
+ // Adding the continuous
+ 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]);
+ }
+ gsl_vector_free(gn);
+}
+
+
+double wgsl_mixed_optim_f(gsl_vector *v, void *params)
+{
+ 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);
+ return mLogLik;
+}
+
+
+/* Compute both f and df together. */
+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);
+}
+
+
+int logistic_mixed_fit(gsl_vector *beta
+ ,gsl_matrix_int *X
+ ,gsl_vector_int *nlev
+ ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc (NULL if not used)
+ ,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;
+
+ //Intializing fix parameters
+ p.X=X;
+ p.Xc=Xc;
+ p.nlev=nlev;
+ p.y=y;
+ p.lambdaL1=lambdaL1;
+ p.lambdaL2=lambdaL2;
+
+ //Initial fit
+ //#ifdef _RPR_DEBUG_
+ mLogLik = wgsl_mixed_optim_f(beta,&p);
+ //fprintf(stderr,"#Initial -log(Lik(0))=%lf\n",mLogLik);
+ //#endif //_RPR_DEBUG
+
+ 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*/
+
+ 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]);
+
+#ifdef _RPR_DEBUG_
+ //mLogLik = wgsl_mixed_optim_f(beta,&p);
+ //fprintf(stderr,"#iter %d, -log(Lik(0))=%lf,%lf\n",(int)iter,mLogLik,maxchange);
+#endif //_RPR_DEBUG
+
+ if(maxchange<1E-4)
+ break;
+ }
+
+#ifdef _RPR_DEBUG_
+ //for (int i = 0; i < npar; i++)
+ // fprintf(stderr,"#par_%d= %lf\n",i,beta->data[i]);
+#endif //_RPR_DEBUG
+
+ //Final fit
+ mLogLik = wgsl_mixed_optim_f(beta,&p);
+ //fprintf(stderr,"#Final %d) -log(Lik(0))=%lf, maxchange %lf\n",iter,mLogLik,maxchange);
+
+ gsl_vector_free (tau);
+ gsl_vector_free (stBeta);
+ gsl_vector_free (myG);
+ gsl_matrix_free (myH);
+
+ return 0;
+}
+
+/***************/
+/* Categorical */
+/***************/
+
+// I need to bundle all the data that goes to the function to optimze together.
+typedef struct{
+ gsl_matrix_int *X;
+ gsl_vector_int *nlev;
+ gsl_vector *y;
+ double lambdaL1;
+ 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)
+{
+ int n = y->size;
+ // int k = X->size2;
+ int npar = beta->size;
+ double total = 0;
+ double aux = 0;
+ /* 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);
+ /* /\*#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;
+ /* #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;
+ }
+ 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 number 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));
+ }
+}
+
+
+/* The gradient of f, df = (df/dx, df/dy). */
+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;
+ /* 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;
+ }
+ // total += y->data[i]*Xbetai-log(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;
+ }
+ }
+}
+
+
+/* The Hessian of f */
+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;
+ int npar = beta->size;
+ // Intitialize Hessian out necessary ???
+ 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
+ // 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?
+ }
+ // total += y->data[i]*Xbetai-log(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?
+ }
+ 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?
+ }
+ }
+}
+
+
+double wgsl_cat_optim_f(gsl_vector *v, void *params)
+{
+ 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);
+ return mLogLik;
+}
+
+
+/* Compute both f and df together. */
+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;
+ fix_parm_cat_T p;
+ int npar = beta->size;
+ int iter=0;
+ double maxchange=0;
+
+ //Intializing fix parameters
+ p.X=X;
+ p.nlev=nlev;
+ p.y=y;
+ p.lambdaL1=lambdaL1;
+ p.lambdaL2=lambdaL2;
+
+ //Initial fit
+ //#ifdef _RPR_DEBUG_
+ mLogLik = wgsl_cat_optim_f(beta,&p);
+ //fprintf(stderr,"#Initial -log(Lik(0))=%lf\n",mLogLik);
+ //#endif //_RPR_DEBUG
+
+ 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*/
+
+ 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]);
+
+#ifdef _RPR_DEBUG_
+ mLogLik = wgsl_cat_optim_f(beta,&p);
+ //fprintf(stderr,"#iter %d, -log(Lik(0))=%lf,%lf\n",(int)iter,mLogLik,maxchange);
+#endif //_RPR_DEBUG
+
+ if(maxchange<1E-4)
+ break;
+ }
+
+#ifdef _RPR_DEBUG_
+ //for (int i = 0; i < npar; i++)
+ // fprintf(stderr,"#par_%d= %lf\n",i,beta->data[i]);
+#endif //_RPR_DEBUG
+
+ //Final fit
+ mLogLik = wgsl_cat_optim_f(beta,&p);
+ //fprintf(stderr,"#Final %d) -log(Lik(0))=%lf, maxchange %lf\n",iter,mLogLik,maxchange);
+
+ gsl_vector_free (tau);
+ gsl_vector_free (stBeta);
+ gsl_vector_free (myG);
+ gsl_matrix_free (myH);
+
+ return 0;
+}
+
+
+/***************/
+/* Continuous */
+/***************/
+
+// 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
+ gsl_vector *y;
+ double lambdaL1;
+ double lambdaL2;
+}fix_parm_cont_T;
+
+
+double fLogit_cont(gsl_vector *beta
+ ,gsl_matrix *Xc
+ ,gsl_vector *y
+ ,double lambdaL1
+ ,double lambdaL2)
+{
+ int n = y->size;
+ int npar = beta->size;
+ double total = 0;
+ double aux = 0;
+ /* 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);
+ /* /\*#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;
+ /* #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));
+ }
+ 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));
+ }
+}
+
+
+/* The gradient of f, df = (df/dx, df/dy). */
+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;
+ /* 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 < Kc; ++k)
+ Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm++];
+
+ pn= -( p->y->data[i] - 1/(1 + gsl_sf_exp(-Xbetai)) );
+
+ 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;
+ }
+ }
+}
+
+
+/* The Hessian of f */
+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;
+ int npar = beta->size;
+ gsl_vector *gn = gsl_vector_alloc(npar); /* gn */
+ // Intitialize Hessian out necessary ???
+ 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
+ // 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++];
+
+ pn= 1/(1 + gsl_sf_exp(-Xbetai));
+ // Add a protection for pn very close to 0 or 1?
+ 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);
+ }
+
+ 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;
+ fix_parm_cont_T *p = (fix_parm_cont_T *)params;
+ 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)
+{
+ *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)
+{
+
+ double mLogLik=0;
+ fix_parm_cont_T p;
+ int npar = beta->size;
+ int iter=0;
+ double maxchange=0;
+
+ //Intializing fix parameters
+ p.Xc=Xc;
+ p.y=y;
+ p.lambdaL1=lambdaL1;
+ p.lambdaL2=lambdaL2;
+
+ //Initial fit
+ //#ifdef _RPR_DEBUG_
+ mLogLik = wgsl_cont_optim_f(beta,&p);
+ //fprintf(stderr,"#Initial -log(Lik(0))=%lf\n",mLogLik);
+ //#endif //_RPR_DEBUG
+
+ 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*/
+
+ 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]);
+
+#ifdef _RPR_DEBUG_
+ mLogLik = wgsl_cont_optim_f(beta,&p);
+ //fprintf(stderr,"#iter %d, -log(Lik(0))=%lf,%lf\n",(int)iter,mLogLik,maxchange);
+#endif //_RPR_DEBUG
+
+ if(maxchange<1E-4)
+ break;
+ }
+
+#ifdef _RPR_DEBUG_
+ //for (int i = 0; i < npar; i++)
+ // fprintf(stderr,"#par_%d= %lf\n",i,beta->data[i]);
+#endif //_RPR_DEBUG
+
+ //Final fit
+ mLogLik = wgsl_cont_optim_f(beta,&p);
+ //fprintf(stderr,"#Final %d) -log(Lik(0))=%lf, maxchange %lf\n",iter,mLogLik,maxchange);
+
+ gsl_vector_free (tau);
+ gsl_vector_free (stBeta);
+ gsl_vector_free (myG);
+ gsl_matrix_free (myH);
+
+ return 0;
+}
+
|