/* Efficient local estimation for time-varying deterministic dynamic models with
application to HIV-1 dynamics
filename: VDM_main.c (J.Chen)
*/

#include "stdio.h"
#include "stdlib.h"
#include <math.h>
#define TWOSTAGE 1
#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
#define MAX(i,j) ( (i) > (j) ? (i):(j) )

main()
{ double *data,*x,*y,**z,*xgrid;
double **esti,**esti1,**VE,**w,**mhat;
double xgridmin,xgridmax;
double t,h,hopt,hopt1,hopt2,c0;
float tempp;
int n,p,q,ngrid,Q,m,N,n0,patient,transfer;
int i,j,k, order;
FILE *reg, *in;
double *vector(), **matrix();
void free_matrix(), free_vector(), exit(), derivate_wle1_data(), constband();
char filename[50], filename0[]="VDM_data.txt";
printf("Local estimation for time-varyinf coefficient HIV dynamic model \n");
printf("Which patient (1 -- 5) ?");
scanf("%d", &patient);
printf("Enter OUTPUT filename:");
scanf("%s", filename);

N=98; transfer = 1;
Q=2; p=1; q=1;

if(patient == 1)
{ c0 = 2.09; n=20; n0 = 0; ngrid = n; /* patient 1 */
}
if(patient == 2)
{ c0 = 2.37; n=19; n0 = 20; ngrid = n; /* patient 2 */
}
if(patient == 3)
{ c0 = 2.4; n=20; n0 = 39; ngrid = n; /* patient 3 */
}
if(patient == 4)
{ c0 = 1.69; n=20; n0 = 59; ngrid = n; /* patient 4 */
}
if(patient == 5)
{ c0 = 2.34; n=18; n0 = 79; ngrid = n; /* patient 5 */
}
xgrid=vector(0,ngrid);
x=vector(0,n);
y=vector(0,n);
z=matrix(0,3,0,n);

VE=matrix(0,3,0,ngrid);
esti = matrix(0,3,0,ngrid);
esti1 = matrix(0,3,0,ngrid);
w = matrix(0,N,0,5);
mhat=matrix(0,ngrid,0,3);

/* INITIALIZATION */
in = fopen(filename0, "r");
if(in == NULL)
{ printf("File not open.\n");
exit(0);
} /*real data*/

reg = fopen(filename,"w");

for(i=0; i < N; i++)
{ for(j=0; j < 5; j++)
{ fscanf(in, "%lf", &w[i][j]);
}
}
fclose(in);

for(i=0; i < n; i++)
{ if(transfer == 1) y[i] = log10(w[i+n0][2]);
else y[i] = w[i+n0][2];
x[i] = w[i+n0][1];
z[0][i] = 1.0;
z[1][i] = w[i+n0][3];
}

printf("\n");
printf("Real AIDS data analysis for the patient %d \n", patient);
for(i=0; i < n;i++)
{ printf("i=%d, x=%6.1f, y=%10.4f z0=%6.1f, z1=%6.1f\n",
i,x[i],y[i],z[0][i],z[1][i]);
}
printf("\n");
printf("The estimaed results for the patient %d \n", patient);

for(i=0; i < ngrid; i++) xgrid[i] = x[i];

order = 1;
constband(x, y, n, order, xgrid, ngrid, mhat, &hopt1);
order = 2;
constband(x, y, n, order, xgrid, ngrid, mhat, &hopt2);


derivate_wle1_data(x,y,z,n,p,q,Q,xgrid,ngrid,VE,esti,esti1,hopt1,hopt2,c0);

fprintf(reg, "Local Two-stage estimation for HIV dynamic model\n");
fprintf(reg, "Real AIDS data analysis\n");
fprintf(reg, "The estimaed results for the patient %d \n", patient);
fprintf(reg, "Estimated bandwidths: hopt1=%10.6f; hopt2=%10.6f c=%6.2f \n", hopt1,hopt2,c0);
fprintf(reg, " i x[i] y[i] EV[i] EV'[i] a1(m) a2(m) a1(L) a2(L) a1(Q) a2(Q)\n");
for(i=0; i < n; i++)
{ fprintf(reg,"%d %8.2f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f \n",
i,x[i],y[i],VE[0][i],VE[1][i],esti[0][i],esti1[0][i],esti[1][i],esti1[1][i],esti[2][i],esti1[2][i]);
}

}

 

void derivate_wle1_data(x,y,z,n,p,q,Q,xgrid,ngrid,VE,esti,esti1,hopt1,hopt2,c)
double *x, *y,**z,*xgrid,**VE,**esti,**esti1,hopt1,hopt2,c;
int n,ngrid,p,q,Q;
#define E 0.00001
#define C 1.71877
{ int i, j, k, m, k1,m1, r, nt,P, P0, order, now0,varying;
double **cof1,**l1,**ll1,*a1,**c1;
double **cof2,**l2,**ll2,*a2,**c2;
double **VE1,**VE2,*yV,**u;
double x0,tmp,tmp0,tmp1,tmp2,tmp3,kh,h0,h1;
double su0,su1,su2,su3,rss,rss0,h,hopt;
double *vector(),**matrix();
void free_vector(),free_matrix(),gaussj(),constband();

P=p+1+q+1;
P0 = P;
cof1 = matrix(0,Q+1,0,ngrid);
l1 = matrix(0,Q+1,0,2);
ll1 = matrix(0,Q+1,0,Q+1);
a1 = vector(0,Q+1);
c1 = matrix(0,Q+1,0,Q+1);

cof2 = matrix(0,P0,0,ngrid);
l2 = matrix(0,P0,0,2);
ll2 = matrix(0,P0,0,P0);
a2 = vector(0,P0);
c2 = matrix(0,P0,0,P0);

VE1 = matrix(0,2,0,n);
VE2 = matrix(0,2,0,n);
yV = vector(0,n);
u=matrix(0,ngrid,0,3);

varying = 1;
su0 = 0.0;
su1 = 0.0;
for(i=1; i < n; i++)
{ su0 += x[i];
su1 += x[i] * x[i];
}
su0 = su0/ (double)(n);
su1 = (su1 - (double)(n) * su0 * su0 )/ (double)(n-1);
su3 = (double)(n);
su2 = pow(su3,-0.2);
h0= 2.34 * sqrt(su1) * su2;
su2 = pow(su3,-0.333);
h1 = sqrt(su1) * su2;
/* One-step: The locally quadrate estimator by fitting Y(t) = V(t) + e */

for(r = 1; r< Q+1; r++)
{ if(r==1) h=0.7*hopt1;
if(r==2) h=0.7*hopt2;
for(j=0; j < n; j++)
{ x0 = x[j];
for(k=0; k< r+1; k++)
{ a1[k] = 0.0;
for(m=0; m < r+1; m++) c1[k][m] = 0.0;
}
for(i=0; i < n; i++)
{ tmp0 = x[i] - x0;
tmp = (x[i] - x0)/h;
if(fabs(tmp) < 1.0)
{ kh = 3.0*(1.0 - tmp*tmp)/(4.0*h);
for(k=0; k < r+1; k++)
{ a1[k] += kh* y[i] * pow(tmp0,(double)k);
for(m=k; m < r+1; m++)
{ c1[k][m] += kh*pow(tmp0,(double)k)*pow(tmp0,(double)m);
}
}
}
}

for(k = 0; k < r; k++)
{ for(m = k+1; m < r+1; m++)
{ tmp = c1[k][m];
c1[m][k] = tmp;
}
}
for(k=0; k < r+1; k++)
{ l1[k][0] = a1[k];
l1[k][1] = 0.0;
for(m=0; m < r+1; m++)
{ ll1[k][m] = c1[k][m];
}
}
gaussj(ll1,r+1,l1,2);

for(k=0; k< r+1; k++) cof1[k][j] = l1[k][0];

if(r==1) /* local linear estimators V(t_j)=VE[0][j] and V'(t_j)=VE[1][j]*/
{ VE1[0][j] = cof1[0][j];
VE1[1][j] = cof1[1][j];
}
if(r==2) /* local quadratic estimators V(t_j)=VE2[0][j] and V'(t_j)=VE2[1][j]*/
{ VE2[0][j] = cof1[0][j];
VE2[1][j] = cof1[1][j];
}
}
}

/* Two-step estimator: Given estimators {V(t_i),V'(t_i)}={VE[0][i],VE[1][i]} */
/* 2.2, Fitting a time-varying coefficients */
if(varying == 1 )
{ h = hopt1;
P = p+1+q+1;

for(r=0; r< 3; r++)
{ if(r==0) /* r=0 mean local mixed estimator u(t) */
{ for(i=0; i < n; i++)
{ VE[0][i] = VE1[0][i];
VE[1][i] = VE2[1][i];
}
}
if(r==1) /* r=1 mean local linear imputed estimator u_1(t) */
{ for(i=0; i < n; i++)
{ VE[0][i] = VE1[0][i];
VE[1][i] = VE1[1][i];
}
}
if(r==2) /* r=2 mean local linear imputed estimator u_2(t) */
{ for(i=0; i < n; i++)
{ VE[0][i] = VE2[0][i];
VE[1][i] = VE2[1][i];
}
}

for(i=0; i < n; i++) yV[i] = VE[1][i] + c * VE[0][i];

order = 1;
constband(x, yV, n, order, xgrid, ngrid, u,&hopt);
hopt2 = hopt;
h = hopt;

for(j=0; j < ngrid; j++)
{ x0 = xgrid[j];
for(k=0; k < P; k++)
{ a2[k] = 0.0;
for(m=0; m < P; m++) c2[k][m] = 0.0;
}
for(i=0; i < n; i++)
{ tmp0 = x[i] - x0;
tmp = (x[i] - x0)/h;
if(fabs(tmp) < 1.0)
{ kh = 3.0*(1.0 - tmp*tmp)/(4.0*h);
for(k=0; k < p+1; k++) a2[k] += kh*yV[i]*z[0][i]*pow(tmp0,(double)k);
for(k=p+1; k < p+1+q+1; k++)
{ k1 = k-(p+1);
a2[k] += kh*yV[i]*z[1][i]*pow(tmp0,(double)k1);
}

for(k=0; k < p+1; k++)
{ for(m=k; m < p+1; m++)
{ c2[k][m] += kh*z[0][i]*z[0][i]*pow(tmp0,(double)k)*pow(tmp0,(double)m);
}
for(m=p+1; m < p+1+q+1; m++)
{ m1=m-(p+1);
c2[k][m] += kh*z[0][i]*z[1][i]*pow(tmp0,(double)k)*pow(tmp0,(double)m1);
}

}
for(k=p+1; k < p+1+q+1; k++)
{ for(m=k; m < p+1+q+1; m++)
{ k1 = k-(p+1); m1 = m-(p+1);
c2[k][m] += kh*z[1][i]*z[1][i]*pow(tmp0,(double)k1)*pow(tmp0,(double)m1);
}
}
}
}

c2[0][0] += 1.0/(5.0*h);
c2[1][1] += 0.2*h/5.0;
c2[2][2] += 1.0/(5.0*h);
c2[3][3] += 0.2*h/5.0;

for(k = 0; k < P-1; k++)
{ for(m = k+1; m < P; m++)
{ tmp = c2[k][m];
c2[m][k] = tmp;
}
}

for(k=0; k < P; k++)
{ l2[k][0] = a2[k];
l2[k][1] = 0.0;
for(m=0; m < P; m++) ll2[k][m] = c2[k][m];
}

gaussj(ll2,P,l2,2);
cof2[0][j] = l2[0][0];
cof2[1][j] = l2[p+1][0];
}
if(r==0)
{ for(j=0; j < n; j++)
{ esti[0][j] = cof2[0][j];
esti1[0][j] = cof2[1][j];
}
}
if(r==1)
{ for(j=0; j < n; j++)
{ esti[1][j] = cof2[0][j];
esti1[1][j] = cof2[1][j];
}
}
if(r==2)
{ for(j=0; j < n; j++)
{ esti[2][j] = cof2[0][j];
esti1[2][j] = cof2[1][j];
}
}
} /* end of time-varying */
printf(" i x[i] V[i] EV[i] EV'[i] a1(m) a2(m); a1(L) a2(L) a1(Q) a2(Q)\n");
for(i=0; i < n; i++)
{ printf("%d %8.3f %10.4f %10.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f \n",
i,x[i],y[i],VE1[0][i],VE1[1][i],esti[0][i],esti1[0][i],esti[1][i],esti1[1][i],esti[2][i],esti1[2][i]);
}

}

for(i=0; i < n; i++)
{ VE[0][i] = VE1[0][i]; VE[1][i] = VE1[1][i];
}
free_matrix(cof1,0,Q+1,0,ngrid);
free_matrix(l1,0,Q+1,0,2);
free_matrix(ll1,0,Q+1,0,Q+1);
free_vector(a1,0,Q+1);
free_matrix(c1,0,Q+1,0,Q+1);

free_matrix(cof2,0,P0,0,ngrid);
free_matrix(l2,0,P0,0,2);
free_matrix(ll2,0,P0,0,P0);
free_vector(a2,0,P0);
free_matrix(c2,0,P0,0,P0);

free_matrix(VE1,0,2,0,ngrid);
free_matrix(VE2,0,2,0,ngrid);
free_vector(yV,0,n);
u=matrix(u,0,ngrid,0,3);
}

 

void constband(x, y, n, order, xgrid, ngrid, mhat, hopt)
double *x, *y, *xgrid, **mhat, *hopt;
int n,order, ngrid;
#define LARGE 1.0e17
#define SPAN 1.1 /*maximum number of iteration*/
#define IUP 3 /*checking conv of iteration*/
#define MIN(i,j) ( (i) > (j) ? (j):(i) )
#define MAX(i,j) ( (i) > (j) ? (i):(j) )
#define ABS(i) ( (i) >= 0 ? (i):-(i))

/* OPTIMAL BANDWIDTH AND REGRESSION ESTIMATION
Input : x :n-column vector containing the covariates;
y : n-column vector with the responses;
order : 0 = local constant approximation
1 = local linear approximation,
2 = local quadratic approximation,
3 = local cubic approximation,
etc. .....,
xgrid : ngrid-column vector with the grid points in
which the final regression estimator mhat is
calculated;
Output : mhat estimated regression curve, bias and var at grid points;
hopt: the optimal bandwidth.
*/

{ FILE *reg;
int i, ii, iup, ismooth, j, nrow, nm, a;
double **beta, **betatmp, *SE, *RSC;
double RSS, h, hmin, hmax,hopt1,MASE, MSEnew, MSEold, xx, xmin, xmax, tmp, sigma;
double hopt2, tmp1, MASEopt, *vector(), **matrix();
void free_vector(), free_matrix();
void sort(), lls();
/*MEMORY ALLOCATION */
nm = order+8;
beta=matrix(0,nm,0,1);
betatmp = matrix(0,nm, 0, ngrid);
SE=vector(0,nm); RSC = vector(0,nm);

xmin=x[0]; xmax=x[0];
for(i=0; i < n; i++)
{if(x[i] < xmin) xmin = x[i];
if(x[i] > xmax) xmax = x[i];
}
hmin = (xmax - xmin)*(order+2.0)/n; /* exp(-0.8*log((double) n)); */
hmax = (xmax - xmin)/2.0;

/* a = AORDER; approximation order*/

if(TWOSTAGE == 1) /* TWO-STAGE PLUGIN */
{ a = AORDER; h = (2.0 +NU)* hmin;}
else
{ a = 0; h = hmin;}

hopt1 = h; iup = 0;
MSEnew = LARGE; MASEopt = LARGE;
/*SEARCH FOR OPTIMAL BANDWIDTH FOR LOCAL order+a FIT */
while(iup < IUP && h < hmax)
{ MASE = 0.0;
for(ii = 0; ii < ngrid; ii++) /* for each grid point */
{xx = xgrid[ii];
/* nrow = # of effective data points */
for(nrow = 0, i=0; i < n; i++) if(fabs(xx - x[i]) < h) nrow++;
if(nrow > 2*(order+a)) /* NOT Nearly singular case */
{ lls(x,y,n,order+a,xx,h, beta, SE, &sigma, RSC, 1);
MASE += RSC[0];
}
else{iup=0; goto lab1;}
}
MSEold = MSEnew;
MSEnew = MASE;
if(MSEnew > MSEold) iup++; else iup = 0;
if(MSEnew < MASEopt) {MASEopt = MSEnew; hopt1 = h;}

lab1: h *= SPAN;
} /*end-loop through h*/
hopt1 *= ADJUST;
hopt2 = hopt1;

if( TWOSTAGE == 1)
{
h = (1.0 + NU)*hmin; hopt2 = h;
iup = 0; MSEnew = LARGE; MASEopt = LARGE;

/* ************* ESTIMATE CURVATIVE FIRST ********** */
for(ii =0; ii < ngrid; ii++)
{ xx = xgrid[ii];
lls(x,y,n,order+a,xx,hopt1, beta, SE, &sigma, RSC, 0);
/*(order+a)-fit */
for(i=0; i <= order+a; i++) betatmp[i][ii] = beta[i][0];
}
/*END priliminary fit to learn beta_2 so that it can be passed */
while(iup < IUP && h < hmax)
{MASE=0.0;
for(ii=0; ii < ngrid; ii++)
{ xx = xgrid[ii];
for(i=0; i <= order+a; i++) beta[i][1] = betatmp[i][ii];
/*END priliminary fit to learn beta_2 so that it can be passed
to compute the bias NOW TRY LOCAL order-FIT,
nrow = # of effective data points */
for(nrow = 0, i=0; i < n; i++) if(fabs(xx - x[i]) < h) nrow++;
if(nrow <= order+a) /* Nearly singular case */
{iup=0; goto lab2;}
else
{ lls(x,y,n,order,xx,h, beta, SE, &sigma, RSC, a);
tmp1 = sigma*sigma;
tmp = tmp1*tmp1;
MASE += (beta[NU][1]*beta[NU][1] + SE[NU]*SE[NU])/tmp;
}
} /*end loop through ii */

MSEold = MSEnew;
MSEnew = MASE;
if(MSEnew > MSEold) iup++; else iup = 0;
if(MSEnew < MASEopt) {MASEopt = MSEnew; hopt2 = h;}

lab2: h *= SPAN;
} /*end loop through h */
/* hopt2 = h * exp(-(IUP+1)*log(SPAN)); */
/* optimal bandwidth for local order-fit*/
} /* end-if with TWOSTAGE == 1 */
/* NOW COMPUTE COMPUTE THE ESTIMATED CURVE */


for(ii=0; ii < ngrid; ii++)
{ xx = xgrid[ii];
lls(x,y,n,order,xx,hopt2, beta, SE, &sigma, RSC, a);
mhat[ii][0] = beta[NU][0];
mhat[ii][1] = beta[NU][1];
mhat[ii][2] = SE[NU];
} /*end the loop through ii */

hopt[0] = hopt2;
free_matrix(beta,0,nm, 0, 1);
free_matrix(betatmp,0,nm,0,ngrid);
free_vector(SE,0,nm);
free_vector(RSC, 0, nm);
}

 

/* ************* 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);
}
 

/* 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);
}
 

/* *******************************LLS.C ************************ */
void lls(x,y,n,p,xx,h, beta, SE, sigma, RSC, a)
double *x, *y,xx,h, **beta, *SE,*sigma, *RSC;
int n, p, a;
/*Do locally weighted least-squared by using local polynomial of order p
Input: x is n-vector of design
y is n-vector of the response
xx is the point where the estimated value to be returned
h is the bandwidth and the kernel is Epnechnikov;
n and p explained above
beta: (p+a) by 2 matrix.
beta[.][1] estimated coefficients to be used to estimate
bias: (p+a)-vector;
a: approximation order in the plub-in procedure and
> 0, return regression coefficients along with others
= 0, return only estimated regression coefficient.

Output: beta[.][0]: estimated regression coefficients (0 -- p)
beta[.][1]: estimated bias according to (2.7) (0 -- p).
SE: p-column vector of estimated SE-values of
the estimated regression coefficients for the
first response vector.
t = beta/SE
sigma: the estimated standard error of linear model.
RSC: Extrended cross-validated values computed as follows.
*/
{
int i,j,l,nrow, nmax;
double RSS, tmp, tmp1, *Xdes,*Y, **H, **Ha, *w ;
double sum, sumw, trace;
double **matrix(), *vector(), *s, *S;
void free_matrix(), free_vector();

/*MEMORY ALLOCATION*/
nmax = MAX(p+a, 2*p+2); /* max moments to be computed */
Xdes = vector(0,n); Y = vector(0,n);
H=matrix(0,p,0,p); Ha = matrix(0,p,0,p);
w = vector(0, n); s = vector(0,nmax);
S = vector(0, nmax);

/* DETERMINING THE LOCAL DATA POINTS */
nrow=0; /* Effective number of data points*/
sumw =0.0; /* Total weight */
for(i=0; i < n; i++)
{tmp = (xx - x[i])/h; tmp *=tmp;
if(tmp < 1.0)
{w[nrow] = 1.0 - tmp; /* Epanechnikov kernel */
sumw += w[nrow];
Xdes[nrow] = x[i] - xx;
Y[nrow] = y[i];
nrow++;
}
}

/* COMPUTE s[j] = \sum w[i] (X[i]-x)^j,
S[j] = \sum w[i]^2 (X[i]-x)^j, beta[j][0] = X'Wy[j] */

for(j=0; j <= nmax; j++) {s[j]=0.0; S[j]=0.0;}
for(j=0; j <= p; j++) beta[j][0] =0.0;
for(i=0; i < nrow; i++)
{ tmp = 1.0; /*tmp = Xdes[i]^j */
for(j=0; j <= nmax; j++)
{tmp1 = w[i]*tmp;
s[j] += tmp1;
S[j] += w[i]*tmp1;
if(j <= p) beta[j][0] += tmp1*Y[i];
tmp *= Xdes[i];
}
}

/*computing H= X'WX, Ha = X'W^2X*/
for(i=0; i <=p; i++)
for(j=0; j <=p; j++)
{ H[i][j] = s[i+j]; Ha[i][j] = S[i+j]; }
/* compute beta[.][1] according to (2.7) */
if(a > 0)
{
for(i=0; i <= p; i++)
{ tmp = 0.0;
for(j=p+1; j <= p+a; j++)
if(j+i <= p+a) tmp += beta[j][1] * s[j+i];
beta[i][1] = tmp;
}
}
gaussj(H,p+1,beta,2); /* solve H b = beta;
output: H <- H^{-1}, beta <- H^{-1} beta */

if(a > 0)
{
for(RSS=.0, i=0; i < nrow; i++) /* compute weighted residual sum of squares
RSS = sum_1^n (y_i - \hat{y}_i)^2 w_i
for the 0th column of Y */
{
sum=0.0; tmp=1.0; /* tmp = Xdes[i]^l */
for(l=0; l <= p; l++)
{ sum += tmp * beta[l][0] ; tmp *= Xdes[i];}
RSS += (Y[i] - sum) * (Y[i] - sum)*w[i];
}
trace = 0.0; /*trace = tr((X'WX)^{-1} X'W^2X )*/
for(i=0; i <= p ; i++)
{ tmp = 0.0;
for(j=0; j<i; j++) tmp += H[i][j] * Ha[j][i];
trace += H[i][i]*Ha[i][i] + 2.0 * tmp;
}
trace = sumw - trace; /*trace = tr(W) - tr((X'WX)^{-1} X'W^2X )*/
RSS = RSS / trace; /* standardarized residual sum of squares */

 

for (l = 0; l <= p; l++)
{ sum = 0.0; /*lth element of (X'WX)^{-1}(X'W^2X)(X'WX)^{-1}*/
for(i=0; i <= p; i++)
{
tmp = 0.0;
for(j=0; j<i; j++) tmp += Ha[i][j] * H[j][l];
sum += H[i][l]*(Ha[i][i]*H[i][l]+2.0 * tmp);
}
SE[l] = sqrt(RSS *sum);
RSC[l] = sum;
}

RSC[0] = RSS * (1.0 + (p+1)*RSC[0]);

sigma[0] = sqrt(RSS);
} /* end-if */

free_vector(Xdes, 0,n); free_vector(Y, 0,n);
free_matrix(H,0,p,0,p); free_matrix(Ha,0,p,0,p);
free_vector(w, 0, n); free_vector(s, 0,nmax);
free_vector(S, 0, nmax);

}