/* Modeling Continuous Auxiliary Covariate Data in the
   Generalized Linear Mixed Models Using the Kernel Smoother
   filename: code.c.  (6/8/2009) By Chen, J       */
/* &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
         Simulated Examples
         
   1,  Logistic regression:  
        P(Y_{ij}=1|X_{ij},u_{i})= e^{alpha+beta*X_{ij}+u_{i} }/(1+e^{alpha+beta*X_{ij}+u_i}
                                                                    j=1,n_j; i=1,n
       Where:  X --- normal distribution N(0,1),
               u_i -- random effects and follows Unif(a,b)
               Z = g(X)+ e  is surrogate variable
               e --- Unif(-sigm_k,sigm_k) distribution   
                                                                                                                                              
   2, Simulation examples         
                                                                                                                                                  
      Input: 1) n  ---- Total sampling size
                v  ---- Total validation sampling size
                s  ---- Number of the groups
                        For each group i, n1_i --- validation number and
                                          n2_i --- novalidation number
                Nsim ---- Number of simulation

             2) alpha ---- coefficient
                beta  ---- coefficient
                sigm_T --- variance of random effect
                sigm_k --- interval  value of errors
                                                                              
   3, Fuctions: The program can run simutancely 4 method, That are
      1),  MLE for Validation data;
      2),  New's MLE for missing and mismeasure data;
      3),  MLE for Validation data without random effects;
      4),  Method for missing without random effects.

  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
 */

#include <stdio.h>
#include <math.h>
#include "Print.c"


#define TWOSTAGE 1                           /* 1 means yes, otherwise no */
#define CONTROL 1
#define AORDER 2
#define NU     0
#define ADJUST (NU == 0? (TWOSTAGE == 1 ? .7776:.8941):(TWOSTAGE == 0 ? .7639:.7643) )
#define NORMIAL         1.96               /* confidence level */
main()
{
   double alpha, beta, sigmat, sigmae, rho0, tmp, tmp1, tmp2, tmp0;
   double *y,**x, **x0, *z,*u,*u0,*mu,*e, *yy, **w, *z0, **mx; 
   double **coef, **mcoef, **vcoef, *rho, **mse, **mv, **var, *sigm, *cov95;
   double su0, su1, su2, su3, h, kh,*xe;
   int n, s, v, p, ps, Nsim, Nsim0,sk, sk1, sk0,skv,P, control, control0;
   int nsim, i, j, k, l, md, md1, md2, md3, md4, md5, first;
   int idum, Idum[3], iexample;
   int *ni, *niv, *cnum, *cnum0, *ivector();
   FILE *reg, *in;
   double *vector(), **matrix(), ran1(), gasdev();
   void sort(), mlem(), free_matrix(), free_vector(), exit();
   char filename[15], filename0[15];

   ps = 15; 
   P = 20; 
   printf("Comparing some methods in semiparameter model with missing data for logistic regression \n"); 
   printf("Which Example (0 -- 1) ?\n Example 0 is real data example, Example 1 is simulation:");
   scanf("%d", &iexample);
 
   if(iexample == 0)
     { printf("Enter sample size n, group number s, covariate number p, P, sigmat \n");
       scanf("%d %d %d %d %lf", &n, &s, &p, &P, &sigmat);
       printf("Whether print data set (yes=1) \n");
       scanf("%d", &control);
       printf("Enter input filename:");
       scanf("%s", filename0);
       printf("Enter OUTPUT filename: ");
       scanf("%s", filename);
       Nsim = 1;         
     }
   else
     {
       printf("Enter sample size n, vali-size v, group number s ( s=>1 ), and Nsim \n");
       scanf("%d %d %d %d", &n, &v, &s, &Nsim);
       printf("Enter parameters: alpha, beta, sigm_T, sigm_k, surrogate (1-3) \n");
       scanf("%lf %lf %lf %lf %d", &alpha, &beta, &sigmat, &sigmae, &first);  
       printf("Enter OUTPUT filename:");
       scanf("%s", filename);
       reg = fopen(filename, "w"); 
       p = 1;  
     }
 
   printf("There are 4 methods in the progrem, you can choose  some of them \n");
   printf("Enter 1, means all of methods; otherwords choose methods(1--4):\n");
   scanf("%d", &md);
   
   if(md != 1) 
    { printf("Method 1 is MLE for validation data (Yes=1) ");
      scanf("%d", &md1);
      printf("Method 2 is New's MLE for missing or mismeasure data (Yes=1)");
      scanf("%d", &md2);
      printf("Method 3 is naive estimator by replaced true as surrogate varibale (Yes=1)");
      scanf("%d", &md3);
      printf("Method 4 is Pepe's method for missing data without random effects (Yes=1)");
      scanf("%d", &md4);
    }

   if(md == 1)
     { md1 = 1;
       md2 = 1;
       md3 = 1;
       md4 = 1;
     }       

    /* MEMORY ALLOCATION */
 ps = p + s + 2; 
 y=vector(0,n);
 yy=vector(0,n);
 x=matrix(0,n,0,p);
 x0=matrix(0,n,0,p);
 mx=matrix(0,p,0,2);
 xe=vector(0,n);

 u=vector(0,s);
 u0=vector(0,s);
 z=vector(0,n);
 z0=vector(0,n);
 e=vector(0,n);

 ni = ivector(0,s+1);
 niv = ivector(0,s+1);
 rho = vector(0,s+1);
 sigm = vector(0,s);

 coef = matrix(0,ps,0,5);
 mcoef = matrix(0,ps,0,5);
 vcoef = matrix(0,ps,0,5);
 mse = matrix(0,ps,0,5);
 var = matrix(0,ps,0,5);
 mv = matrix(0,ps,0,5);
 mu = vector(0,s);
 sigm = vector(0,s+1);

 cov95 = vector(0,5);
 cnum = ivector(0,5);
 cnum0 = ivector(0,5);
 w = matrix(0,n,0,P);

 switch(iexample)
     {
      case 0: in = fopen(filename0, "r");       
         reg = fopen(filename,"w");
         for(i=1; i < s+1; i++)
           { fscanf(in, "%d %d ",ni+i, niv+i);
           }
       
         for(i=0; i < n; i++)
           { for(j=0; j < P; j++)
               { fscanf(in, "%lf", &w[i][j]); 
               }
           }
 
         for(i=0; i < n; i++)  
           { y[i] = w[i][2];
             yy[i] = w[i][2]; 
             z[i] = w[i][10];
             x[i][0] = w[i][6];
             x[i][1] = w[i][8];  
             x[i][2] = w[i][9]-1.0;
             x[i][3] = w[i][10];
             x[i][4] = w[i][14];
             x[i][5] = w[i][15];
             x[i][6] = w[i][16];
             x[i][7] = w[i][17];
             x[i][8] = w[i][18];
             x[i][9] = w[i][12];
           }

    for(k=0; k<p; k++)
       { mx[k][0] = 0.0;
         mx[k][1] = 0.0;
       } 
    for(k=0; k<p; k++) 
       { for(i=0; i<n; i++)
           { mx[k][0] += x[i][k];
             mx[k][1] += x[i][k] * x[i][k];
           }
       }
    for(k=0; k<p; k++)
       { mx[k][0] = mx[k][0]/(double)n;
         tmp =  mx[k][1]/(double)n - mx[k][0] * mx[k][0];
         mx[k][1] = sqrt(tmp);
       } 

    idum = -10; 
    for(k=0; k < s; k++)
       { u[k] = sigmat * gasdev(&idum);
       }

           
    su0 = 0.0;
    su1 = 0.0;
    for(i=0; i< n; i++)
      { su0 += z[i];
        su1 += z[i] * z[i];
      }
    su0 = su0/(double)(n);
    su1 = (su1-(double)(n) * su0 * su0 ) / (double)(n-1);
    su3 = (double)(n);
    su2 = pow(su3,-0.3333);
    h = 2.0 * sqrt(su1) * su2;

    sk0 = 0;
    for(i=0; i < s; i++)
       { sk0 += ni[i];
         sk = sk0 + ni[i+1];   
         skv = sk0 + niv[i+1];
         for(j= sk0; j < skv; j++)
           { su0 = 0.0; su1 = 0.0; 
             for(l= sk0; l < skv; l++)
               { tmp = (z[l]-z[j])/h;
                 if(fabs(tmp) <1.0)
                    { kh = 3.0 * (1.0 - tmp*tmp)/(4.0*h);
                      su0 += kh;
                      su1 += kh * x[l][0];
                    }   
               }
             xe[j] = su1/su0;
           }
       }
   
    for(i=0; i < n; i++)
      { printf("i=%d  %7.4f %7.4f %7.4f %7.4f \n",i,h,xe[i],x[i][0],z[i]);
      }              
 

    if(control == 1)
       { fprintf(reg, " i, center,  y[i],  z[i],  x[i][0], x[i][1] --- \n");
         for(i=0; i < n; i++)
           {  fprintf(reg,"%d %5.1f %7.2f %7.2f", i, w[i][1], y[i], z[i]);
              for(j=0; j < p; j++)
                {  fprintf(reg,"%7.2f  ", x[i][j]);
                }
              fprintf(reg,"%10.2f", xe[i]);  
              fprintf(reg, "\n");
           }
       }
 
         printf("Real data analysis \n");

         fprintf(reg, "Likelihood estimations with missing covariate data for a mixed effects model ession\n");
         fprintf(reg, "Real data: n=%d, s=%d, p=%d ps=%d P=%d sigma=%8.3f \n", n, s, p, ps,P,sigmat);
      
         mlem(y,yy,x,z,n,s,p,ps,ni,niv,coef,var,md1,md2,md3,md4,md5,u);
       
         for(i=0; i<s; i++)  rho[i]  = (double)(niv[i]) / (double)(ni[i]);
 
         for(j=0; j < 4; j++)       
           { for(i=0; i < p+1; i++) 
               { var[i][j] = sqrt( var[i][j] );    
               }
           }
         fprintf(reg, "Descriptive data:\n"); 
         fprintf(reg, "Factors(i)  MEAN      S.D \n");   
         for(i=0; i< p; i++)
            { fprintf(reg, "%d %16.12f %16.12f \n", i,mx[i][0],mx[i][1]);
            }

         fprintf(reg, "Estimate results\n");
     
        if(md1 == 1)
         { fprintf(reg, "1, MLE for validation data :\n");
           fprintf(reg, "i   Estimation of coef(i)   SE(coef(i)) \n");   
           for(i=0; i< p+1; i++)  
             { fprintf(reg, "%d %16.12f %16.12f \n", i, coef[i][0],var[i][0]);
             }
         }
  
        if(md2 == 1)
         { fprintf(reg, "2, New Method:\n");
           fprintf(reg, "i   Estimation of coef(i)   SE(coef(i)) \n");
           for(i=0; i< p+1; i++)
             { fprintf(reg, "%d   %16.12f   %16.12f \n", i, coef[i][1], var[i][1]);
             }           
         }
 
        if(md3 == 1)
         { fprintf(reg, "3, MLE for valitation without considering random effects:\n");
           fprintf(reg, "i   Estimation of coef(i)   SE(coef(i)) \n");
           for(i=0; i< p+1; i++)
             { fprintf(reg, "%d   %16.12f   %16.12f \n", i, coef[i][2],var[i][2]);
             }
         }
     
        if(md4 == 1)
         { fprintf(reg, "4, CW' method without considering random effects:\n");
           fprintf(reg, "i   Estimation of coef(i)   SE(coef(i)) \n");
           for(i=0; i< p+1; i++)
             { fprintf(reg, "%d   %16.12f   %16.12f \n", i,coef[i][3],var[i][3]);
             }      
         }

        fprintf(reg, "\n");
        fprintf(reg," i   n[i]   nv[i]  u[i]   u[i]-Mle   u[i]-New \n");  
        for(i=1; i < s+1; i++)
        fprintf(reg, "%5d  %5d %5d  %16.12f   %16.12f  %16.12f\n",i,ni[i],niv[i], u[i+1], coef[p+i][0],coef[p+i][1]);
    

          exit(1);
     
       case 1:                  
            for(i=0; i < n; i++)
              { tmp = 1;
                   y[i]= alpha;
              }
            break;

	   default:  printf("Out of Example list:\n");
           exit(1);
      }


ni[0] = 0; niv[0] = 0; 
for(k=1; k < s+1; k++) ni[k] = n/s;
for(k=1; k < s+1; k++) niv[k] = v/s;
for(k=0; k < s; k++) 
  { idum = -10*k*k-2+k;
    sigm[k] = sigmae; 
    /* sigm[k] = ran1(&idum) * sigmae; */
  } 
        /* sigm[k] = sigmae +  ran1(&idum) * (1.0-sigmae); */
 

fprintf(reg, "Likelihood estimations of GLMM with continue surrogate information (code.c) \n");
fprintf(reg, "EXAMPLE %d,  Covariate number p= %d, Total estimated parameters =%d \n",iexample, p, ps);
fprintf(reg, "Sample size n = %d, group number s =%d,  validation size v = %d,  number of simulations = %d \n",
              n, s, v, Nsim);
fprintf(reg, "True parameters: alpha= %8.5f, beta=%8.5f, sigmat=%8.5f, sigmae=%8.5f, auxiary=%d\n",
               alpha, beta, sigmat, sigmae, first);
fprintf(reg,"\n\n");
/*
fprintf(reg, "nsim;  alpha1  beta1; alpha2  beta2; alpha3 beta3; alpha4  beta4\n");
 */

/* *************************************************************
                    SIMULATION STARTS
   ************************************************************  */

 printf("Comparing some methods for GLMM with continue surrogate information\n");
 printf("Simulation Experience \n");
 printf("Example %d:\n",iexample);
                      
 for(j=0; j < 5; j++)
   { for(i=0; i < ps; i++)
      { mcoef[i][j] = 0.0;
        vcoef[i][j] = 0.0;
        mse[i][j] = 0.0;
        mv[i][j] = 0.0;
      }  
   }

 for(k=0; k < s; k++) mu[k] = 0.0;
 for(j=0; j < 5; j++) cnum[j] = 0;
 Nsim0 =0; 


 for(nsim=0; nsim < Nsim; nsim++)
   {
     idum = -10*(nsim+1)*(nsim+1)-2+(nsim+1);
  
     if(s == 1)
       { for(k=0; k < s; k++)
           { u[k] = 0.0;
             mu[k] += u[k];
           } 
       }
     else
       { for(k=0; k < s; k++) 
           { if(nsim == nsim) 
               { u0[k] = sigmat * gasdev(&idum);
                 u[k] = u0[k];
               } 
             else
               { u[k] = u0[k];
               } 
             mu[k] += u[k];
           }
       } 

     sk = 0;
     for(k=0; k < s; k++)
       { sk += ni[k];
         sk1 = sk + ni[k+1];
         for(i= sk; i < sk1; i++)
           {  if(first == 1)
                { x[i][0] =  2.0 * ran1(&idum); 
                  e[i] = sigm[k] * (2.0 * ran1(&idum) - 1.0);
                  z[i] = x[i][0] + e[i];
                }
              if(first == 2)
                { x[i][0] =  2.0 * ran1(&idum);
                  e[i] = sigm[k] * (2.0 * ran1(&idum) - 1.0);
                  z[i] = x[i][0] * x[i][0] + e[i];
                }
              if(first == 3)
                { x[i][0] =  2.0 * ran1(&idum);
                  e[i] = sigm[k] * (2.0 * ran1(&idum) - 1.0);                  
                  z[i] = cos((3.14179/2.0) * x[i][0]) + e[i];
                }
           }
       }

     sk = 0; 
     for(k=0; k < s; k++)  
       {  sk += ni[k];
          sk1 = sk + ni[k+1];
          for(i= sk; i < sk1; i++)  
            { tmp = alpha + beta * x[i][0];
              tmp1 = exp(tmp+u[k])/(1.0+exp(tmp+u[k]));
              tmp2 = exp(tmp)/(1.0+exp(tmp));
              tmp0 = ran1(&idum);              
              if( tmp0 < tmp1)  y[i] = 1.0;
              else y[i]=0.0;
              if( tmp0 < tmp2)  yy[i] = 1.0;
              else yy[i]=0.0;
            }
       } 

     su0 = 0.0;
     for(k=0; k < s; k++) su0 += u[k]*u[k];              
     su0 = sqrt( su0/(double)s );  
     

     sk0 = 0;
     tmp0 = 0.0; tmp1 = 0.0;
     for(i=0; i < s; i++)
        { sk0 += ni[i];
          sk1 = sk0 + ni[i+1];
          skv = sk0 + niv[i+1];
          for(j= sk0; j < skv; j++)
            { tmp = alpha + beta * x[j][0];                  
              tmp0 += y[j]*(tmp + u[i]);
              tmp1 += log(1.0+ exp(tmp + u[i]));
            }
         }
     tmp2=0.5*(double)s*log(2.0*3.14*sigmat*sigmat)+0.5*(double)(s*s)*su0*su0/(sigmat*sigmat);         
     su1 = tmp0 - tmp1 -tmp2;   


   for(i=0; i < s; i++) u[i] = 0.0;  
 
   mlem(y,yy,x,z,n,s,p,ps,ni,niv,coef,var,md1,md2,md3,md4,u,sigmat); 
  
   tmp =0.0;
   for(j=0; j < 5; j++)
     { for(i=0; i < 2; i++) 
          { tmp += coef[i][j];
          }
     }
  
   if(tmp < 10.0)     
    { Nsim0 ++; 
      for(j=0; j < 5; j++)
       { for(i=0; i < ps; i++)   
           { mcoef[i][j] += coef[i][j];
             vcoef[i][j] += coef[i][j] * coef[i][j];
             mv[i][j] += var[i][j];
           } 
         mse[0][j] += (coef[0][j] - alpha) * (coef[0][j] - alpha);
         mse[1][j] += (coef[1][j] - beta) * (coef[1][j] - beta);
       }
    

      for(j=0; j < 5; j++) cnum0[j] = 0.0;
      for(j=0; j < 5; j++)
         {  tmp0 =  1.98 * sqrt( var[1][j] );
            tmp1 = coef[1][j] - beta;
            if( fabs(tmp1) < tmp0 ) 
               { cnum[j]++;
                 cnum0[j] = 1.0;
               }
         }

    } 
    /*
     fprintf(reg,"%d  %9.5f %8.5f %9.5f %8.5f %9.5f %8.5f %9.5f %8.5f %d %10.5f %d %8.5f %d %8.5f %d %8.5f \n",
       nsim,coef[0][0],coef[1][0],coef[0][1],coef[1][1],coef[0][2],coef[1][2],coef[0][3],coef[1][3], cnum0[0],var[1][0],cnum0[1],var[1][1],cnum0[2],var[1][2], cnum0[3],var[1][3]);
     */   
  }  
   
   fprintf(reg,"\n");
   Nsim = Nsim0;
   for(j=0; j < 5; j++) cov95[j] = (double)(cnum[j])/ (double)(Nsim); 	
   for(j=0; j < 5; j++)
       { for(i=0; i < ps; i++)
           { mcoef[i][j] = mcoef[i][j] / (double)(Nsim);
             vcoef[i][j] = (vcoef[i][j] -(double)(Nsim)*mcoef[i][j]*mcoef[i][j])/(double)(Nsim-1);
             vcoef[i][j] = sqrt( vcoef[i][j]);
             mv[i][j] =  sqrt( mv[i][j]/(double)(Nsim) ); 
           }
         mse[0][j] = sqrt( mse[0][j]/(double)(Nsim-1) );
         mse[1][j] = sqrt( mse[1][j]/(double)(Nsim-1) );   
       }
   for(k=0; k < s; k++) mu[k] = mu[k]/ (double)(Nsim);

   rho0 = (double)(v) / (double)(n);

 if(md1 == 1)
  {fprintf(reg, "1, MLE for Validation data :\n");
   fprintf(reg,
"    rho; mean((alpha0,beta0)); var((alpha0,beta0)); MSE((alpha0,beta0); mvar(alpha0,beta0); cov(beta0)\n");
   fprintf(reg," %7.4f %9.5f %9.5f %9.5f %9.5f  %9.5f %9.5f %9.5f %9.5f %9.5f \n",
      rho0, mcoef[0][0], mcoef[1][0], vcoef[0][0], vcoef[1][0],mse[0][0],mse[1][0],  mv[0][0],mv[1][0], cov95[0]);
  }

 if(md2 == 1)
  {fprintf(reg,"\n");
   fprintf(reg, "2, New's MLE for missing data :\n");
   fprintf(reg,
"    rho; mean((alpha0,beta0)); var((alpha0,beta0)); MSE((alpha0,beta0); mvar(alpha0,beta0); cov(beta0)\n");
   fprintf(reg," %7.4f %9.5f %9.5f %9.5f %9.5f  %9.5f %9.5f %9.5f %9.5f %9.5f \n",
      rho0, mcoef[0][1], mcoef[1][1], vcoef[0][1], vcoef[1][1], mse[0][1],mse[1][1], mv[0][1],mv[1][1], cov95[1]);
  }  

 if(md3 == 1)
  {fprintf(reg,"\n");
   fprintf(reg, "3, Naive estimator by replaced surrogate as true variable :\n");
   fprintf(reg,
"    rho; mean((alpha0,beta0)); var((alpha0,beta0)); MSE((alpha0,beta0); mvar(alpha0,beta0); cov(beta0)\n");
   fprintf(reg," %7.4f %9.5f %9.5f %9.5f %9.5f  %9.5f %9.5f %9.5f %9.5f %9.5f \n",
      rho0, mcoef[0][2], mcoef[1][2], vcoef[0][2], vcoef[1][2], mse[0][2],mse[1][2], mv[0][2],mv[1][2],cov95[2]);
  }  
 
 if(md4 == 1)
  {fprintf(reg,"\n");
   fprintf(reg, "4, CW's method for missing data without random effects :\n");
   fprintf(reg,
"    rho; mean((alpha0,beta0)); var((alpha0,beta0)); MSE((alpha0,beta0); mvar(alpha0,beta0); cov(beta0)\n");
   fprintf(reg," %7.4f %9.5f %9.5f %9.5f %9.5f  %9.5f %9.5f %9.5f %9.5f %9.5f \n",
      rho0, mcoef[0][3], mcoef[1][3], vcoef[0][3], vcoef[1][3], mse[0][3],mse[1][3], mv[0][3],mv[1][3],cov95[3]);
  }  
 
   fclose(reg);  
 
}       

/************************** MLEM.C*********************************/  
void  mlem(y,yy,x,z,n,s,p,ps,ni,niv,coef,var,md1,md2,md3,md4,u,sigmat)
double *y, *yy, **x, *z, **coef, **var,*u,sigmat;
int n, s, p, ps, *ni, *niv, md1,md2,md3,md4;
#define E  0.0001
#define C  1.71877

{
  int i, j, i0, j0, r, k, l, k0,l0, k1,l1,m, p1, p2, jj, nt, nv;
  int sk0, sk, skv,sk01,sk1,skv1,sk02,sk2,skv2,num;
  double tmp, tmp0, tmp1, tmp2, tmp3, tmu, tmu1, tmu2, su0, su1, su2, su3, su4, su5, su6;
  double sm0,sm1,sm2,sm3,sm4,sm5,sm6,pv, xxi, xxj;
  double h, *hopt, *hopt0, kh, sumk,fun;
  double **cof,*a0,**c0,**a,**c,**id, **idv, *d, **x0, *dp0, **dp1,**v, **v0, **w, *xz;  
  double *vector(), **matrix();
  void free_vector(), free_matrix();
  void gaussj(), mtv();
 
  cof = matrix(0,ps,0,5);          
  a0 = vector(0,ps);
  c0 = matrix(0,ps,0,ps);
  a = matrix(0,ps,0,2);
  c = matrix(0,ps,0,ps);
  id = matrix(0,n,0,n);
  idv = matrix(0,n,0,n);
  d = vector(0,n);
  x0 = matrix(0,n,0,p+1);
  xz = vector(0,n);

  dp0 = vector(0,n);
  dp1 = matrix(0,n,0,ps);

  v = matrix(0,ps,0,ps);
  v0 = matrix(0,ps,0,ps);
  w = matrix(0,n,0,ps);
 
  hopt = vector(0,s);
  hopt0 = vector(0,s);

  p1 = p+1;
  nv = 0;
  for(i=0; i<s; i++)  nv += niv[i];
  pv = (double)(nv)/(double)(n);
  
  for(i=0; i<n; i++)
    { x0[i][0] = 1.0;
      for(j=0; j<p; j++)
        { tmp = x[i][j];
          x0[i][j+1] = tmp;
        }
    }
 
  su0 = 0.0;
  su1 = 0.0;
  for(i=0; i<n; i++)
    { su0 += z[i]; 
      su1 += z[i] * z[i];
    }
  su0 = su0/ (double)(n);
  su1 = (su1 - (double)(n) * su0 * su0 )/ (double)(n-1);
  su3 = (double)(n)/(double)(s);
  su2 = pow(su3,-0.3333);
  h = 2.0 *  sqrt(su1) * su2;
  for(i=0; i < s; i++) hopt[i] = 2.5 * h;
  
  for(i=0; i<n; i++)
    { x0[i][0] = 1.0;
      for(j=0; j<p; j++)
        { tmp = x[i][j];
          x0[i][j+1] = tmp;
        }
    }

  sk0 = 0;
  for(i=0; i < s; i++)
    { sk0 += ni[i];
      sk = sk0 + ni[i+1];    
      skv = sk0 + niv[i+1];
      for(j = sk0; j < skv; j++) xz[j] = x0[j][1];         
      for(j = skv; j < sk; j++) xz[j] = z[j];
    }


       /*  1, MLE Method for validation data and mixed effects model  */

if(md1 == 1)
 { for(i=0; i < p1; i++) cof[i][0] = 0.0;
   cof[p1+s][0] = 0.5;
   for(i=0; i < s; i++)
     { k = p1 + i;
       cof[k][0] = 0.2;
     }

   if(s == 1)
     { p2 = p1;
     }
   else
     { p2 = p1 + s;
     }
  
   for(m=1; m <= 1000; m++)
     { for(i0=0; i0 < p2; i0++)
         { a0[i0] = 0.0;
           for(j0=0; j0 < p2; j0++)
             { c0[i0][j0] = 0.0;
             }
         }                       
       
       if(s == 1)
         { for(i=0; i < s; i++)  cof[p1+i][0] = 0.0;
           tmu1 = 1.0;
         }
       else 
         { /*tmu1 = cof[p1+s][0];*/    
            tmu1 = sigmat;   
         }
       tmu2 = tmu1 * tmu1;

       sk0 = 0;
       for(i=0; i < s; i++)
         { sk0 += ni[i];
           sk = sk0 + ni[i+1];
           skv = sk0 + niv[i+1];
           for(j= sk0; j < skv; j++)    
              { su0 = 0.0;
                for(r=0; r < p1; r++)
                  { su0 += x0[j][r] * cof[r][0];
                  }
                tmp0 = exp( su0 + cof[p1+i][0] + u[i]);
                tmp1 = tmp0/(1.0 + tmp0);
                tmp2 = tmp1/(1.0 + tmp0);
                for(i0=0; i0 < p1; i0++)
                  { su1 = (y[j] - tmp1) * x0[j][i0];
                    a0[i0] += su1;
                    for(j0=i0; j0 < p1; j0++)
                       { su2 = - tmp2 * x0[j][i0] * x0[j][j0];
                         c0[i0][j0] += su2;
                       }
                  }             
              }
         } 
              
       if(s != 1)
         { sk0 = 0;
           for(i=0; i < s; i++)
            { sk0 += ni[i];
              sk = sk0 + ni[i+1];     
              skv = sk0 + niv[i+1];
              k = p1 + i;
              for(j= sk0; j < skv; j++)
                { su0 = 0.0;
                  for(r=0; r < p1; r++)
                    { su0 += x0[j][r] * cof[r][0];
                    }
                  tmp0 = exp( su0 + cof[p1+i][0] + u[i]);
                  tmp1 = tmp0/(1.0 + tmp0);
                  tmp2 = tmp1/(1.0 + tmp0);
                  a0[k] += y[j] - tmp1;
                  for(i0=0; i0 < p1; i0++)
                    { c0[i0][k] += - tmp2*x0[j][i0];
                    }
                  c0[k][k] +=  -tmp2;  
                }
              a0[k] += - cof[k][0]/tmu2;
              c0[k][k] +=  - 1.0/tmu2;
              for(j0=k+1; j0< p2; j0++) c0[k][j0] = 0.0;
            }            
         }                          

       for(i0=0; i0 < p2-1; i0++)
         { for(j0=i0+1; j0 < p2; j0++)
             { tmp = c0[i0][j0];
               c0[j0][i0] = tmp;
             }
         }

       for(i=0; i< p2; i++)
         { a[i][0] = a0[i];
           a[i][1] = 0.0;
           for(j=0; j < p2; j++)
             { c[i][j] = c0[i][j];
             }
         }

       gaussj(c,p2,a,2);
       for(j=0; j < p2; j++)
         { cof[j][0] += -a[j][0];
         }

       
       if(s != 1)
         { su3 = 0.0;  
           for(i=0; i < s; i++)
             { k = p1+i;
               su3 += cof[k][0] * cof[k][0];
             }
           cof[p1+s][0] = sqrt( su3/(double)(s) );
         }
                 /* compute the likelihood function */
       sk0 = 0; 
       tmp0 = 0.0; tmp1 = 0.0; 
       for(i=0; i < s; i++)
         { sk0 += ni[i];   
           sk = sk0 + ni[i+1];
           skv = sk0 + niv[i+1];
           for(j= sk0; j < skv; j++)
              { su0 = 0.0;
                for(r=0; r < p1; r++)
                  { su0 += x0[j][r] * cof[r][0];
                  }  
                tmp0 += y[j]*( su0+ cof[p1+i][0] );
                tmp1 += log(1.0+ exp(su0 + cof[p1+i][0]));
              }
         }
       tmp2 = 0.5* (double)s * log(2.0*3.14*tmu2) + 0.5*(double)(s*s)*cof[p1+s][0]*cof[p1+s][0]/tmu2;                       
       fun = tmp0 - tmp1 -tmp2;

       su4 = 0.0;
       for(j=0; j < p+1; j++)
          { su4 += a[j][0];
          }
       if(fabs(su4) < E) break;
     }

     for(j=0; j < p2; j++) var[j][0] = -c[j][j];
     for(j=0; j < ps; j++) coef[j][0] = cof[j][0];
 }   
                                  

      /*  2, New's Method for missing data and mixed effects model  */
if(md2 == 1)
 { for(i=0; i < p1; i++) cof[i][1] = 0.0;
   cof[p1+s][1] = 0.25;
   for(i=0; i < s; i++)  
     { k = p1 + i;
       cof[k][1] = 0.01;
     }

   if(s == 1)
     { p2 = p1;
     }
   else
     { p2 = p1 + s;
     }

   for(m=1; m <= 100; m++)
     { for(i0=0; i0< p2; i0++) 
         { a0[i0] = 0.0;
           for(j0=0; j0< p2; j0++) 
             { c0[i0][j0] = 0.0;
             }
         }

       if(s == 1)
         { for(i=0; i < s; i++)  cof[p1+i][1] = 0.0;
           tmu1 = 1.0;
         }
       else
         { /*tmu1 = cof[p1+s][1];*/
           tmu1 = sigmat;
         }
       tmu2 = tmu1 * tmu1;

       sk0 = 0; 
       for(i=0; i < s; i++)
         { sk0 += ni[i];
           sk = sk0 + ni[i+1];
           skv = sk0 + niv[i+1];
           for(j= sk0; j < skv; j++)
              { su0 = 0.0;      
                for(r=0; r<p1; r++)            
                  { su0 += x0[j][r] * cof[r][1];
                  }
                tmp0 = exp( su0 + cof[p1+i][1] + u[i]);
                tmp1 = tmp0/(1.0 + tmp0);
                tmp2 = tmp1/(1.0 + tmp0);
                for(i0=0; i0 < p1; i0++)
                  { su1 = (y[j] - tmp1) * x0[j][i0];
                    a0[i0] += su1;   
                    for(j0=i0; j0 < p1; j0++)
                       { su2 = - tmp2*x0[j][i0]*x0[j][j0];
                         c0[i0][j0] += su2;
                       }
                  } 
              }
         }

    if(s != 1) 
     { sk0 = 0; 
       for(i=0; i < s; i++)
          { sk0 += ni[i];
            sk = sk0 + ni[i+1];
            skv = sk0 + niv[i+1];
            k = p1 + i;
            for(j= sk0; j < skv; j++)  
              { su0 = 0.0;
                for(r=0; r < p1; r++)
                  { su0 += x0[j][r] * cof[r][1];
                  }
                tmp0 = exp( su0+cof[p1+i][1] + u[i]);
                tmp1 = tmp0/(1.0 + tmp0);
                tmp2 = tmp1/(1.0 + tmp0);
                a0[k] += y[j] - tmp1;
                for(i0=0; i0 < p1; i0++)
                 {  c0[i0][k] += - tmp2 * x0[j][i0];
                 }
                c0[k][k] +=  - tmp2;
              }
            a0[k] += - cof[k][1]/tmu2;
            c0[k][k] +=  - 1.0/tmu2; 
            for(j0=k+1; j0< p2; j0++) c0[k][j0] = 0.0;
          }
     }
          
             /*  Imputing the missing data  */

       for(i=0; i < n; i++) d[i] = 0.0;
       sk0 = 0;      
       for(i=0; i < s; i++)
          { h = hopt[i];
            sk0 += ni[i];
            sk = sk0 + ni[i+1];
            skv = sk0 + niv[i+1];
            for(j= skv; j < sk; j++)
              { sumk = 0.0;
                num = 0;
                su0 = z[j]; 
                for(l= sk0; l < skv; l++)
                   { tmp = (z[l]-z[j])/h;
                     if(fabs(tmp) <1.0) 
                       { kh = 3.0 * (1.0 - tmp*tmp)/(4.0*h);
                         sumk += kh; 
                       }
                   }

                for(l= sk0; l < skv; l++)
                   { tmp = (z[l]-z[j])/h;
                     if(fabs(tmp) <1.0)
                       { kh = 3.0*(1.0 - tmp*tmp)/(4.0*h);
                         id[j][l] = kh/sumk;
                       }
                     else   
                       { id[j][l] = 0.0;
                       }
                   }                 
             
                for(l= sk0; l < skv; l++)
                   { su0 = 0.0;
                     for(r=0; r < p1; r++)
                       { if( r == 1 )
                           { su0 += x0[l][r] * cof[r][1];
                           }
                         else
                           { su0 += x0[j][r] * cof[r][1];
                           }   
                       }
                     tmp0 = exp( y[j] * (su0+cof[p1+i][1] + u[i]) );
                     tmp1 = 1.0/( 1.0+exp(su0+cof[p1+i][1]) );
                     d[j] += tmp0 * tmp1 * id[j][l];
                   }

                dp0[j] = d[j];
              }
          }
    
       sk0 = 0;  
       for(i=0; i < s; i++) 
         { sk0 += ni[i];  
           sk = sk0 + ni[i+1];   
           skv = sk0 + niv[i+1];   
           for(j= skv; j < sk; j++)
             { for(i0=0; i0 < p1; i0++)
                 { for(j0=i0; j0 < p1; j0++)  
                     { sm0 = 0.0; sm1 = 0.0; sm2 = 0.0;
                       sm3 = 0.0; sm4 = 0.0; sm5 = 0.0;
                       for(l= sk0; l < skv; l++)
                          { su0 = 0.0;
                            for(r=0; r < p1; r++)
                                { if( r == 1 )
                                   { su0 += x0[l][r] * cof[r][1];
                                   }
                                  else
                                   { su0 += x0[j][r] * cof[r][1];
                                   }
                                }
                 
                            if(i0 == 1)
                              { xxi = x0[l][i0];
                              } 
                            else
                              { xxi = x0[j][i0];
                              }
                            if(j0 == 1)
                              { xxj = x0[l][j0];
                              }
                            else
                              { xxj = x0[j][j0];
                              }
   
                            tmp0 = exp( su0+cof[p1+i][1] + u[i]);
                            tmp1 = tmp0/(1.0 + tmp0);
                            tmp2 = 1.0/(1.0 + tmp0);
                            tmp3 = tmp2/(1.0 + tmp0);
                            su2 = exp( y[j]* ( su0+cof[p1+i][1] + u[i]) )* tmp2 * id[j][l];
             
                            sm0 += (y[j] - tmp1) * su2;
                            sm1 += (y[j] - tmp1) * su2 * xxi;
                            sm2 += (y[j] - tmp1) * su2 * xxj;
                            sm4 += su2;
                            if(y[j] == 1.0)
                              { sm3 += (1.0 - tmp0) * tmp3 * su2 * xxi * xxj;
                              }
                            else
                              { sm3 += (tmp0 -1.0) * tmp0 * tmp3 * su2 * xxi * xxj;
                              }
                          }
                       c0[i0][j0] += sm3/sm4-sm1*sm2/(sm4*sm4);
                     }
                   a0[i0] += sm1/sm4;
                   dp1[j][i0] = sm1;
                 }                  
             }
         } 

    if(s != 1)
     { sk0 = 0;
       for(i=0; i < s; i++)   
         { sk0 += ni[i];
           sk = sk0 + ni[i+1];
           skv = sk0 + niv[i+1];
           k0 = p1 + i;
           for(j= skv; j < sk; j++)
             { for(i0=0; i0<p1; i0++)
                 { sm0 = 0.0; sm1 = 0.0; sm2 = 0.0;
                   sm3 = 0.0; sm4 = 0.0; sm5 = 0.0;
                   for(l= sk0; l < skv; l++)
                      { su0 = 0.0;
                        for(r=0; r < p1; r++)
                          { if(r == 1)
                              { su0 += x0[l][r] * cof[r][1];
                              }
                            else
                              { su0 += x0[j][r] * cof[r][1];
                              }
                          }
                            
                        if(i0 == 1)
                          { xxi =  x0[l][i0];
                          } 
                        else
                          { xxi =  x0[j][i0];
                          }

                        tmp0 = exp( su0 + cof[p1+i][1] + u[i]);
                        tmp1 = 1.0/(1.0 + tmp0);
                        tmp2 = tmp1/(1.0 + tmp0);
                        su1 = y[j] - tmp0 * tmp1;
                        su2 = exp( y[j] * ( su0+cof[p1+i][1] + u[i]) ) * tmp1 * id[j][l];
                        sm0 += su1 * su2;
                        sm1 += su1 * su2 * xxi;
                        if(y[j] == 1.0)
                          { sm2 += (1.0 - tmp0)*tmp2*su2;
                            sm3 += (1.0 - tmp0)*tmp2*su2*xxi;
                          }
                        else
                          { sm2 += (tmp0 - 1.0)*tmp0*tmp2*su2;
                            sm3 += (tmp0 - 1.0)*tmp0*tmp2*su2*xxi;
                          }
                      }
   
                    c0[i0][p1+i] += sm3/d[j]- sm0*sm1/(d[j]*d[j]);
                    c0[p1+i][p1+i] += sm2/d[j]-sm0*sm0/(d[j]*d[j]);
                 }
                dp1[j][p1+i] = sm0;
                a0[p1+i] += sm0/d[j];
                for(j0=k0+1; j0< p2; j0++) c0[k0][j0] = 0.0;
             }
         }  
     }

       for(i0=0; i0 < p2-1; i0++)
         { for(j0=i0+1; j0 < p2; j0++)
             { tmp = c0[i0][j0];
               c0[j0][i0] = tmp;
             }
         }

       for(i=0; i<p2; i++)
         { a[i][0] = a0[i];
           a[i][1] = 0.0;  
           for(j=0; j<p2; j++)
             { c[i][j] = c0[i][j];
             }
         }
       gaussj(c,p2,a,2);

       for(j=0; j<p2; j++)
          { cof[j][1] += -a[j][0];
          }

    if(s != 1)
     { su3 = 0.0;
       for(i=0; i < s; i++)
         { k = p1+i;
           su3 += cof[k][1] * cof[k][1];
         }
       cof[p1+s][1] = sqrt( su3/(double)(s) );
     }

       su4 = 0.0;
       for(j=0; j<p+1; j++)
          { su4 += a[j][0];
          }        

/*
      printf("m=%d  h=%7.4f %7.4f %7.4f %7.4f \n",m,h,cof[0][1],cof[1][1],su4);
 */               
 
       if(fabs(su4) < E) break;
        if(fabs(su4) > 10) break; 

     }   
                                                    
   /*  Computing  variance of estimator  */
   
    sk0 = 0.0;
    for(i=0; i < s; i++)
       { h = hopt[i];
         sk0 += ni[i];
         sk = sk0 + ni[i+1];
         skv = sk0 + niv[i+1];
         for(j= sk0; j < skv; j++)
           { num = 0;
             sumk = 0.0;
             for(l= skv; l < sk; l++)
               { tmp = (z[l]-z[j])/h;
                 if(fabs(tmp) < 1.0) 
                    { kh = 3.0*(1.0 - tmp*tmp)/(4.0*h);
                      sumk += kh;
                    }
               } 

             for(l= skv; l < sk; l++)
               { tmp = (z[l]-z[j])/h;
                 if(fabs(tmp) < 1.0)
                    { kh = 3.0*(1.0 - tmp*tmp)/(4.0*h);
                      idv[j][l] = kh / sumk;
                    }
                 else
                    { idv[j][l] = 0.0;
                    }
               }
           }
       }
 
     sk0 = 0;
     for(i=0; i < s; i++)
       { sk0 += ni[i];
         sk = sk0 + ni[i+1];
         skv = sk0 + niv[i+1];
         for(j= sk0; j < skv; j++)
           { for(i0=0; i0 < p1; i0++)
               { sm0 = 0.0; sk01 = 0;
                 for(l= skv; l < sk; l++)
                    { su0 = 0.0;
                      for(r=0; r < p1; r++)
                         { if( r == 1 )
                             { su0 += x0[j][r] * cof[r][1];
                             }
                           else
                             { su0 += x0[l][r] * cof[r][1];
                             }
                         }

                      if(i0 == 1)
                         { xxi =  x0[j][i0];
                         } 
                      else
                         { xxi =  x0[l][i0];
                         }

                      tmp0 = exp( su0+ cof[p1+i][1] + u[i]);
                      tmp1 = tmp0/(1.0 + tmp0);
                      tmp2 = 1.0/(1.0 + tmp0);
                      tmp3 = tmp2/(1.0 + tmp0);
                      su1 = exp( y[l] * ( su0 + cof[p1+i][1] + u[i]) )* tmp2;
                      su2 = (y[l] - tmp1) * su1 * xxi;
                      sm0 += ( su2 - dp1[l][i0]*su1/dp0[l] ) * idv[j][l]/dp0[l]; 
                    }

                 w[j][i0] = sm0;
               }
             if(s != 1)
               { for(i0=p1; i0 < p1+s; i0++)
                   { sm0 = 0.0; sk01 = 0;            
                     k = i0 - p1;
                     sk01 += ni[k]; 
                     sk1 = sk01 + ni[k+1];
                     skv1 = sk01 + niv[k+1];
                     for(l= skv1; l < sk1; l++)
                       { su0 = 0.0;
                         for(r=0; r < p1; r++)
                            { if( r == 1 )
                                { su0 += x0[j][r] * cof[r][1];
                                }
                              else
                                { su0 += x0[l][r] * cof[r][1];
                                }
                            }
           
                         if(i0 == 1)
                            { xxi =  x0[j][i0];
                            } 
                         else
                            { xxi =  x0[l][i0];
                            }
                         tmp0 = exp( su0 + cof[i0][1] + u[i] );
                         tmp1 = tmp0/(1.0 + tmp0);
                         tmp2 = 1.0/(1.0 + tmp0); 
                         tmp3 = tmp2/(1.0 + tmp0);
                         su1 = exp( y[l] * ( su0 + cof[i0][1] + u[i]) )* tmp2;
                         su2 = (y[l] - tmp1) * su1;
                         sm0 += ( su2 - dp1[l][i0]*su1/dp0[l]) * idv[j][l]/dp0[l];
                       }
                     w[j][i0] = sm0;
                   }
               }
           }
       }
 	
    p2 = p1;
    for(i0=0; i0< p2; i0++)
       { for(j0=0; j0< p2; j0++)
           { v[i0][j0] = 0.0;
           }
       } 
    for(i0=0; i0< p2; i0++)
       { for(j0=0; j0< p2; j0++)
           { su1 = 0.0; 
             su2 = 0.0;
             sk0 = 0;
             for(i=0; i < s; i++)
               { sk0 += ni[i];
                 sk = sk0 + ni[i+1];
                 skv = sk0 + niv[i+1];
                 for(j= sk0; j < skv; j++)
                   { su1 += w[j][i0];
                     su2 += w[j][j0];
                   }
               }

             su1 = su1/(double)(nv);
             su2 = su2/(double)(nv);
             sk0 = 0;  
             for(i=0; i < s; i++)
               { sk0 += ni[i];   
                 sk = sk0 + ni[i+1];
                 skv = sk0 + niv[i+1];
                 for(j= sk0; j < skv; j++)
                   { v[i0][j0] += (w[j][i0]-su1) * (w[j][j0]-su2)/(double)(nv-1);
                   }
               }
           }  
       }

    mtv(c,v,v0,p2,p2,p2);  
    mtv(v0,c,v,p2,p2,p2);
 
    for(i=0; i < p2; i++) var[i][1] = - c[i][i] + (1.0-pv)*(1.0-pv) * v[i][i]/pv;
    for(j=0; j < ps; j++) coef[j][1] = cof[j][1];
 }  


     /*  3, MLE  Method for validation data without random effects   */

/*
for(i=0; i < n; i++) y[i] = yy[i]; 

if(md3 == 1)
 { for(i=0; i < p1; i++) cof[i][2] = 0.0;
   for(m=1; m <= 1000; m++)
     { for(i0=0; i0< p1; i0++)
         { a0[i0] = 0.0;
           for(j0=0; j0< p1; j0++)
             { c0[i0][j0] = 0.0;
             }
         }  

       sk0 = 0;
       for(i=0; i < s; i++)
         { sk0 += ni[i];
           sk = sk0 + ni[i+1];
           skv = sk0 + niv[i+1];
           for(j= sk0; j < skv; j++)
              { su0 = 0.0;
                for(r=0; r < p1; r++)
                  { su0 += x0[j][r] * cof[r][2];
                  }
                tmp0 = exp( su0 + u[i] );
                tmp1 = tmp0/(1.0 + tmp0);
                tmp2 = tmp1/(1.0 + tmp0); 
                for(i0=0; i0 < p1; i0++)
                  { su1 = (y[j] - tmp1) * x0[j][i0];
                    a0[i0] += su1;
                    for(j0=i0; j0 < p1; j0++)
                       { su2 = - tmp2*x0[j][i0]*x0[j][j0];
                         c0[i0][j0] += su2;
                       }
                  }
              }
         }
          
       for(i0=0; i0 < p; i0++)
         { for(j0=i0+1; j0 < p1; j0++)  
             { tmp = c0[i0][j0]; 
               c0[j0][i0] = tmp;
             }
         }
                  
       for(i=0; i< p1; i++)
         { a[i][0] = a0[i];
           a[i][1] = 0.0;
           for(j=0; j< p1; j++)
             { c[i][j] = c0[i][j];
             }
         }    
 
       gaussj(c,p1,a,2);
             
       for(j=0; j < p1; j++)
          { cof[j][2] += -a[j][0];
          }

       su4 = 0.0;
       for(j=0; j < p1; j++)
          { su4 += a[j][0];
          }
       if(fabs(su4) < E) break;   
     }
                  
     for(j=0; j < p1; j++) var[j][2] = -c[j][j];
     for(j=0; j < ps; j++) coef[j][2] = cof[j][2];
 }
  */ 
      
      /*  4, Method for missing data */

if(md4 == 1)
 { for(i=0; i < p1; i++) cof[i][3] = 0.0;
   for(m=1; m <= 100; m++)
     { for(i0=0; i0< p1; i0++)
         { a0[i0] = 0.0;   
           for(j0=0; j0< p1; j0++)
             { c0[i0][j0] = 0.0;
             }
         }
               
       sk0 = 0;
       for(i=0; i < s; i++)
         { sk0 += ni[i];
           sk = sk0 + ni[i+1];
           skv = sk0 + niv[i+1];
           for(j= sk0; j < skv; j++)
              { su0 = 0.0;
                for(r=0; r < p1; r++)
                  { su0 += x0[j][r] * cof[r][3];
                  }
                tmp0 = exp(su0 + u[i]);
                tmp1 = tmp0/(1.0 + tmp0);
                tmp2 = tmp1/(1.0 + tmp0);
                for(i0=0; i0 < p1; i0++)
                  { a0[i0] += (y[j] - tmp1) * x0[j][i0];
                    for(j0=i0; j0 < p1; j0++)
                       { c0[i0][j0] += - tmp2 * x0[j][i0] * x0[j][j0];
                       }
                  }
              }  
         }
       
            /*  Imputing the missing data  */

       for(i=0; i < s; i++)
          { hopt0[i] = 1.5 * hopt[i]; 
          }
       sk0 = 0;
       for(i=0; i < s; i++)
          { h = hopt0[i];
            sk0 += ni[i];
            sk = sk0 + ni[i+1];
            skv = sk0 + niv[i+1];
            for(j= skv; j < sk; j++)
              { sumk = 0.0;
                num = 0;
                su0 = z[j];
                for(l= sk0; l < skv; l++)
                   { tmp = (z[l]-z[j])/h;
                     if(fabs(tmp) <1.0)
                       { kh = 3.0*(1.0 - tmp*tmp)/(4.0*h);
                         sumk += kh;
                       }
                   } 
       
                for(l= sk0; l < skv; l++)
                   { tmp = (z[l]-z[j])/h;
                     if(fabs(tmp) <1.0)
                       { kh = 3.0*(1.0 - tmp*tmp)/(4.0*h);
                         id[j][l] = kh/sumk;
                       }
                     else
                       { id[j][l] = 0.0;
                       }
                   }
              }
          }

       sk0 = 0;
       for(i=0; i < s; i++)
         { sk0 += ni[i];
           sk = sk0 + ni[i+1];
           skv = sk0 + niv[i+1];
           for(j= skv; j < sk; j++)
             { for(i0=0; i0 < p1; i0++)
                 { for(j0=i0; j0 < p1; j0++)
                     { sm0 = 0.0; sm1 = 0.0; sm2 = 0.0;
                       sm3 = 0.0; sm4 = 0.0; sm5 = 0.0;
                       for(l= sk0; l < skv; l++)
                          { su0 = 0.0;
                            for(r=0; r < p1; r++)
                                { if( r == 1 )
                                   { su0 += x0[l][r] * cof[r][3];
                                   }
                                  else
                                   { su0 += x0[j][r] * cof[r][3];
                                   }
                                }
                            if(i0 == 1)
                              { xxi =  x0[l][i0];
                              } 
                            else
                              { xxi =  x0[j][i0];
                              }

                            if(j0 == 1)
                              { xxj =  x0[l][j0];
                              } 
                            else
                              { xxj =  x0[j][j0];
                              }

                            tmp0 = exp( su0+u[i] );
                            tmp1 = tmp0/(1.0 + tmp0);
                            tmp2 = 1.0/(1.0 + tmp0);
                            tmp3 = tmp2/(1.0 + tmp0);
                            su2 = exp( y[j] * ( su0 + u[i]) )* tmp2 * id[j][l];
                             
                            sm0 += (y[j] - tmp1) * su2;
                            sm1 += (y[j] - tmp1) * su2 * xxi;
                            sm2 += (y[j] - tmp1) * su2 * xxj;
                            sm4 += su2;
                            if(y[j] == 1.0)
                               { sm3 += (1.0 - tmp0) * tmp3 * su2 * xxi * xxj;
                               }
                            else
                               { sm3 += (tmp0 - 1.0) * tmp0 * tmp3 * su2 * xxi * xxj;
                               }
                          }
                         
                       c0[i0][j0] += sm3/sm4-(sm1*sm2)/(sm4*sm4);
                     }
                   a0[i0] += sm1/sm4;
                   dp0[j] = sm4;
                   dp1[j][i0] = sm1;
                 }
             }
         }

       for(i0=0; i0 < p; i0++)
         { for(j0=i0+1; j0 < p1; j0++)
             { tmp = c0[i0][j0];
               c0[j0][i0] = tmp;
             }
         }
                         
       for(i=0; i < p1; i++)
         { a[i][0] = a0[i];   
           a[i][1] = 0.0;
           for(j=0; j < p1; j++)
             { c[i][j] = c0[i][j];
             }
         }
       gaussj(c,p1,a,2);
                    
       for(j=0; j < p1; j++)
          { cof[j][3] += -a[j][0];
          }
               
       su4 = 0.0;
       for(j=0; j < p1; j++)
          { su4 += a[j][0];  
          }

       if(fabs(su4) < E) break;
        if(fabs(su4) >10) break; 
     }

  /*  Comput variance of estimator  */

    sk0 = 0.0;                  
    for(i=0; i < s; i++)
       { h = hopt0[i];
         sk0 += ni[i];
         sk = sk0 + ni[i+1];
         skv = sk0 + niv[i+1];
         for(j= sk0; j < skv; j++)
           { num = 0;
             sumk = 0.0;
             for(l= skv; l < sk; l++)
               { tmp = (z[l]-z[j])/h;
                 if(fabs(tmp) < 1.0)
                    { kh = 3.0*(1.0 - tmp*tmp)/(4.0*h);
                      sumk += kh;
                    } 
               }
   
             for(l= skv; l < sk; l++)
               { tmp = (z[l]-z[j])/h;
                 if(fabs(tmp) < 1.0)
                    { kh = 3.0*(1.0 - tmp*tmp)/(4.0*h);
                      idv[j][l] = kh / sumk;
                    }
                 else
                    { idv[j][l] = 0.0;
                    }
               }
           }
       }
   
     sk0 = 0;
     for(i=0; i < s; i++)
       { sk0 += ni[i];
         sk = sk0 + ni[i+1];
         skv = sk0 + niv[i+1];
         for(j= sk0; j < skv; j++)
           { for(i0=0; i0 < p1; i0++)
               { sm0 = 0.0; sk01 = 0;
                 for(l= skv; l < sk; l++)
                    { su0 = 0.0;
                      for(r=0; r < p1; r++)
                        { if( r == 1 )
                            { su0 += x0[j][r] * cof[r][3];
                            }
                          else
                            { su0 += x0[l][r] * cof[r][3];
                            }
                        }
                      if(i0 == 1)
                         { xxi =  x0[j][i0];
                         } 
                      else
                         { xxi =  x0[l][i0];
                         }

                      tmp0 = exp( su0 + u[i] );
                      tmp1 = tmp0/(1.0 + tmp0);
                      tmp2 = 1.0/(1.0 + tmp0);
                      tmp3 = tmp2/(1.0 + tmp0);
                      su1 = exp( y[l] * ( su0+u[i] ) )* tmp2;
                      su2 = (y[l] - tmp1) * su1 * xxi;
                      sm0 += ( su2 - dp1[l][i0]*su1/dp0[l] ) * idv[j][l]/dp0[l];
                    }
                 w[j][i0] = sm0;
               }
           }
       }
                     
    for(i0=0; i0< p1; i0++)
       { for(j0=0; j0< p1; j0++)
           { v[i0][j0] = 0.0;
           }
       }
    for(i0=0; i0< p1; i0++)
       { for(j0=0; j0< p1; j0++)
           { su1 = 0.0;
             su2 = 0.0;
             sk0 = 0;
             for(i=0; i < s; i++)
               { sk0 += ni[i];
                 sk = sk0 + ni[i+1];  
                 skv = sk0 + niv[i+1];
                 for(j= sk0; j < skv; j++)
                   { su1 += w[j][i0];
                     su2 += w[j][j0];
                   }
               }
             su1 = su1/(double)(nv);
             su2 = su2/(double)(nv);
             sk0 = 0;
             for(i=0; i < s; i++)
               { sk0 += ni[i];
                 sk = sk0 + ni[i+1];
                 skv = sk0 + niv[i+1];
                 for(j= sk0; j < skv; j++)
                   { v[i0][j0] += (w[j][i0]-su1) * (w[j][j0]-su2)/(double)(nv-1);
                   }
               }
           }
       }
                 
    mtv(c,v,v0,p1,p1,p1);
    mtv(v0,c,v,p1,p1,p1);
    for(i=0; i < p1; i++) var[i][3] = -c[i][i] + (1.0-pv)*(1.0-pv) * v[i][i]/pv;
    for(j=0; j < ps; j++) coef[j][3] = cof[j][3];
 }


  /* 3. Naive estimator */

if(md3 == 1)
 { for(i=0; i < p1; i++) cof[i][2] = 0.0;
   cof[p1+s][2] = 0.25;
   for(i=0; i < s; i++)
     { k = p1 + i;
       cof[k][2] = 0.2;
     }
          
   if(s == 1)
     { p2 = p1;
     }
   else
     { p2 = p1 + s;
     }
           
   for(m=1; m <= 1000; m++)
     { for(i0=0; i0 < p2; i0++)
         { a0[i0] = 0.0;
           for(j0=0; j0 < p2; j0++)
             { c0[i0][j0] = 0.0;
             }
         }
                
       if(s == 1)
         { for(i=0; i < s; i++)  cof[p1+i][2] = 0.0;
           tmu1 = 1.0;
         }
       else
         { tmu1 = cof[p1+s][2];
         }
       tmu2 = tmu1 * tmu1;
               
       sk0 = 0;
       for(i=0; i < s; i++)
         { sk0 += ni[i];
           sk = sk0 + ni[i+1];
           skv = sk0 + niv[i+1];
           for(j= sk0; j < sk; j++)
              { 
                x0[j][1] = xz[j];
                su0 = 0.0;
                for(r=0; r < p1; r++)
                  { su0 += x0[j][r] * cof[r][2];
                  }
                tmp0 = exp( su0 + cof[p1+i][2] + u[i]);
                tmp1 = tmp0/(1.0 + tmp0);
                tmp2 = tmp1/(1.0 + tmp0);
                for(i0=0; i0 < p1; i0++)
                  { su1 = (y[j] - tmp1) * x0[j][i0];
                    a0[i0] += su1;
                    for(j0=i0; j0 < p1; j0++)
                       { su2 = - tmp2 * x0[j][i0] * x0[j][j0];
                         c0[i0][j0] += su2;
                       }
                  }
              }
         }

   if(s != 1)
         { sk0 = 0;
           for(i=0; i < s; i++)
            { sk0 += ni[i];
              sk = sk0 + ni[i+1];
              skv = sk0 + niv[i+1];
              k = p1 + i;
              for(j= sk0; j < sk; j++)
                { x0[j][1] = xz[j];
                  su0 = 0.0;
                  for(r=0; r < p1; r++)
                    { su0 += x0[j][r] * cof[r][2];
                    }
                  tmp0 = exp( su0 + cof[p1+i][2] + u[i]);
                  tmp1 = tmp0/(1.0 + tmp0);
                  tmp2 = tmp1/(1.0 + tmp0);
                  a0[k] += y[j] - tmp1;
                  for(i0=0; i0 < p1; i0++)
                    { c0[i0][k] += - tmp2*x0[j][i0];
                    }
                  c0[k][k] +=  - tmp2;
                }
              a0[k] += - cof[k][2]/tmu2;
              c0[k][k] +=  - 1.0/tmu2;
              for(j0=k+1; j0< p2; j0++) c0[k][j0] = 0.0;
            }
         }

     for(i0=0; i0 < p2-1; i0++)
         { for(j0=i0+1; j0 < p2; j0++)
             { tmp = c0[i0][j0];
               c0[j0][i0] = tmp;
             }
         }
                    
       for(i=0; i< p2; i++)
         { a[i][0] = a0[i];
           a[i][1] = 0.0;
           for(j=0; j < p2; j++)
             { c[i][j] = c0[i][j];
             }
         }
                     
       gaussj(c,p2,a,2);
       for(j=0; j < p2; j++)
         { cof[j][2] += -a[j][0];
         }
       if(s != 1)
         { su3 = 0.0;
           for(i=0; i < s; i++)  
             { k = p1+i;
               su3 += cof[k][2] * cof[k][2];
             }
           cof[p1+s][2] = sqrt( su3/(double)(s) );
         }
                    
       su4 = 0.0;
       for(j=0; j < p+1; j++)
          { su4 += a[j][0];
          }
       if(fabs(su4) < E) break;   
     }
          
     for(j=0; j < p2; j++) var[j][2] = -c[j][j];
     for(j=0; j < ps; j++) coef[j][2] = cof[j][2];
 }
         


  free_matrix(cof,0,ps,0,5);  
  free_vector(a0,0,ps);
  free_matrix(c0,0,ps,0,ps);
  free_matrix(a,0,ps,0,2);
  free_matrix(c,0,ps,0,ps);

  free_matrix(id,0,n,0,n);
  free_matrix(idv,0,n,0,n); 
  free_vector(d,0,n);
  free_matrix(x0,0,n,0,p+1);
  free_vector(xz,0,n);

  free_matrix(v,0,ps,0,ps);
  free_matrix(v0,0,ps,0,ps);
  free_matrix(w,0,n,0,ps);
  free_vector(dp0,0,n);
  free_matrix(dp1,0,n,0,ps);

  free_vector(hopt,0,s);
  free_vector(hopt0,0,s);
 }


/* *************  MTV.C ***********************************
   Function:     Matrix time vector
      Input:  a: n by k input matrix
              b: k by m input matrix
     Output:  c: n by m matrix,   i.e c = a * b  
  *********************************************************/
void mtv(a,b,c,n,k,m)
double **a,**b,**c;
int n,k,m;
                  
{ int i,j,l;
  double sum;
  void free_vector(), free_matrix();
           
  for(i=0; i < n; i++)
    { for(j=0; j < m; j++)
        { sum = 0.0;
          for(l=0; l < k; l++) sum += a[i][l] * b[l][j];
          c[i][j] = sum;
        }
    }
  
}

/* *************  GAUSSJ.C ************** */
/* #include <math.h>
#include <stdio.h>  */
#define SWAP(a,b) {double temp=(a);(a)=(b);(b)=temp;}
#define TINY 1.0e-12;
void gaussj(a,n,b,m)
     /* Linear equation solution by Gauss-Jordan elimination
        Input: a: n by n input matrix
               b: n by m input matrix containing the m right-hand sides
                  vectors of the equation-system
        In the present version, b is a column vector (m=1). The program
        has to be changed slightly for the general case (m > 1). In this
        case, b should be declared as a matrix and e.g. line 58 should
        look like:
            for (l=1; l<=m; l++) SWAP(b[irow-1][l-1], b[icol-1][l-1])
        Output: a: the input matrix a is replaced be its matrix inverse
                b: the input matrix is replaced by the corresponding set
                   of solution vectors. */

double **a,**b;
int n,m;

{
  int *indxc,*indxr,*ipiv;
  int i,icol,irow,j,k,l,ll, *ivector();
  double big,dum,pivinv;
  void nrerror(), free_ivector();

  indxc=ivector(1,n);
  indxr=ivector(1,n);
  ipiv=ivector(1,n);


  for (j=1; j <= n; j++) ipiv[j]=0;
  for (i=1; i <= n; i++) {
        big=0.0;
        for (j=1; j <= n; j++)
            if (ipiv[j] !=1)
                 for (k=1; k<=n; k++) {
                   if (ipiv[k] == 0) {
                       if (fabs(a[j-1][k-1]) >= big) {
                           big=fabs(a[j-1][k-1]);
                           irow=j;
                           icol=k;
                       }
                   } 
     /*else if (ipiv[k] > 1) nrerror("GAUSSJ: Singular Matrix-1");
      */
      }
        ++(ipiv[icol]);

        if (irow !=icol) {
            for (l=1; l<=n; l++) SWAP(a[irow-1][l-1], a[icol-1][l-1])
            for (l=1; l<=m; l++) SWAP(b[irow-1][l-1], b[icol-1][l-1])
        }
        indxr[i]=irow;
        indxc[i]=icol;
        if (fabs(a[icol-1][icol-1]) == 0.0) a[icol-1][icol-1] = TINY;
        pivinv=1.0/a[icol-1][icol-1];
        a[icol-1][icol-1]=1.0;
        for(l=1;l<=n;l++) a[icol-1][l-1] *=pivinv;
        for(l=1;l<=m;l++) b[icol-1][l-1] *=pivinv;
        for (ll=1;ll<=n;ll++)
            if (ll !=icol) {
                dum=a[ll-1][icol-1];
                a[ll-1][icol-1]=0.0;
                for (l=1;l<=n;l++) a[ll-1][l-1] -= a[icol-1][l-1]*dum;
                for (l=1;l<=m;l++) b[ll-1][l-1] -= b[icol-1][l-1]*dum;
            }
  }
  for (l=n; l>=1;l--) {
      if (indxr[l] !=indxc[l])
          for (k=1;k<=n;k++)
                  SWAP(a[k-1][indxr[l]-1],a[k-1][indxc[l]-1]);
  }
  free_ivector(ipiv,1,n);
  free_ivector(indxr,1,n);
  free_ivector(indxc,1,n);
}

/*       ********************* SORT.C ***************** */
void sort(ra,n)
int n;
double *ra;
/*index of ra is 0:n-1*/
{
        int l,j,ir,i;
        double rra;

        l=(n >> 1)+1;
        ir=n;
        for (;;) {
                if (l > 1)
                        rra=ra[--l-1];
                else {
                        rra=ra[ir-1];
                        ra[ir-1]=ra[0];
                        if (--ir == 1) {
                                ra[0]=rra;
                                return;
                        }
                }
                i=l;
                j=l << 1;
                while (j <= ir) {
                        if (j < ir && ra[j-1] < ra[j]) ++j;
                        if (rra < ra[j-1]) {
                                ra[i-1]=ra[j-1];
                                j += (i=j);
                        }
                        else j=ir+1;
                }
                ra[i-1]=rra;
        }
}

/* ******************** UNIFORM RANDOM NUMBER GENERATOR: RAN1.c ******* */
#include <stdlib.h>
#define M1 259200
#define IA1 7141
#define IC1 54773
#define RM1 (1.0/M1)
#define M2 134456
#define IA2 8121
#define IC2 28411
#define RM2 (1.0/M2)
#define M3 243000
#define IA3 4561
#define IC3 51349

double ran1(idum)
int *idum;

                    /* returns a uniform random deviate between 0.0 and 1.0.
                       Set idum to any negative value to initialize or
                       reinitialize the sequence  */

{
    static long ix1,ix2,ix3;
    static double r[98];
    double temp;
    static int iff=0;
    int j;
    void nrerror();

    if (*idum < 0 || iff == 0) {
         iff=1;
         ix1=(IC1-(*idum)) % M1;
         ix1=(IA1*ix1+IC1) % M1;
         ix2=ix1 % M2;
         ix1=(IA1*ix1+IC1) % M1;
         ix3=ix1 % M3;
         for (j=1;j<=97;j++) {
               ix1=(IA1*ix1+IC1) % M1;
               ix2=(IA2*ix2+IC2) % M2;
               r[j]=(ix1+ix2*RM2)*RM1;
           }
           *idum=1;
    }
    ix1=(IA1*ix1+IC1) % M1;
    ix2=(IA2*ix2+IC2) % M2;
    ix3=(IA3*ix3+IC3) % M3;
    j=1 + ((97*ix3)/M3);
    if (j > 97 || j < 1) nrerror("RAN1: This cannot happen.");
    temp=r[j];
    r[j]=(ix1+ix2*RM2)*RM1;
    return temp;
}

/* ***************** NORMAL RANDOM GENERATOR: GASDEV.C *********** */
double gasdev(idum)
int *idum;
                   /* returns a normally distributed deviate with mean 0
                      and variance 1, using ran1(idum) as the source of uniform
                      deviates    */

{
    static int iset=0;
    static double gset;
    double fac,r,v1,v2;
    double ran1();

    if (iset == 0)  {
        do {
            v1=2.0*ran1(idum)-1.0;
            v2=2.0*ran1(idum)-1.0;
            r=v1*v1+v2*v2;
        } while (r >= 1.0 || r == 0.0);
        fac= sqrt(-2.0*log(r)/r);
        gset=v1*fac;
        iset=1;
        return v2*fac;
    }  else {
         iset=0;
         return gset;
    }
}

/*                    UTILITY PROGRAMS  */
/* ************************************************************************* */
/*        *************** IVECTOR.C  **********    */
int *ivector(nl, nh)
int nl, nh;
{  int *v;
   v = (int *)malloc((unsigned) (nh-nl+1) * sizeof(int));
   return v-nl;
}

/*        *************** VECTOR.C  **********    */
double *vector(nl, nh)
int nl, nh;
{  double *v;
   v = (double *)malloc((unsigned) (nh-nl+1) * sizeof(double));
   return v-nl;
}


/*        ***************MATRIX.C  **********    */
double **matrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;

        /* Allocates a double matrix with range [nrl..nrh][ncl..nch]    */

{
    int i;
    double **m;

    m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
    m -=nrl;

    for(i=nrl;i<=nrh;i++) {
         m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
         m[i] -= ncl;
    }
    return m;
}

/*        ***************FREE_MATRIX.C  **********    */
void free_matrix(m,nrl,nrh,ncl,nch)
double **m;
int nrl,nrh,ncl,nch;
                     /* Frees a matrix allocated with matrix */

{
    int i;

    for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
    free((char*) (m+nrl));
}

/*        ***************FREE_IVECTOR.C  **********    */
void free_ivector(v,nl,nh)
int *v,nl,nh;
                       /* Frees an int vector allocated by ivector() */

{
     free((char*) (v+nl));
}

/*        ***************FREE_VECTOR.C  **********    */
void free_vector(v,nl,nh)
double *v;
int nl,nh;
                       /* Frees a double vector allocated by vector() */

{
     free((char*) (v+nl));
}


/*    *********************NRERROR.C ******************** */
void nrerror(error_text)
     /* prints an error message */

char error_text[];
{
   fprintf(stderr,"Numerical recipes run-time error...\n");
   fprintf(stderr,"%s\n",error_text);
   fprintf(stderr,"...now exiting to system...\n");
   exit(1);
}