1 /*
2  *   Copyright (c) 1996-2001 Lucent Technologies.
3  *   See README file for details.
4  *
5  *
6  *   Post-fitting functions to compute the local variance and
7  *   influence functions. Also the local degrees of freedom
8  *   calculations for adaptive smoothing.
9  */
10 
11 #include "local.h"
12 
13 extern double robscale;
14 
15 /*
16   vmat() computes (after the local fit..) the matrix
17   M2  = X^T W^2 V X.
18   M12 = (X^T W V X)^{-1} M2
19   Also, for convenience, tr[0] = sum(wi) tr[1] = sum(wi^2).
20 */
vmat(lf,des,M12,M2,tr)21 void vmat(lf, des, M12, M2, tr)
22 lfit *lf;
23 design *des;
24 double *M12, *M2, *tr;
25 { INT i, p, nk, ok;
26   double link[LLEN], h, ww;
27   p = des->p;
28   setzero(M2,p*p);
29 
30   nk = -1;
31 
32   /* for density estimation, use integral rather than
33      sum form, if W^2 is programmed...
34   */
35   if ((lf->mi[MTG]<=THAZ) && (lf->mi[MLINK]==LLOG))
36   { switch(lf->mi[MKER])
37     { case WGAUS: nk = WGAUS; h = des->h/SQRT2; break;
38       case WRECT: nk = WRECT; h = des->h; break;
39       case WEPAN: nk = WBISQ; h = des->h; break;
40       case WBISQ: nk = WQUQU; h = des->h; break;
41       case WTCUB: nk = W6CUB; h = des->h; break;
42       case WEXPL: nk = WEXPL; h = des->h/2; break;
43     }
44   }
45 
46   tr[0] = tr[1] = 0.0;
47   if (nk != -1)
48   { ok = lf->mi[MKER]; lf->mi[MKER] = nk;
49 /* compute M2 using integration. Use M12 as work matrix. */
50     (des->itype)(des->xev, M2, M12, lf, des->cf, h);
51     lf->mi[MKER] = ok;
52     if (lf->mi[MTG]==TDEN) multmatscal(M2,lf->dp[DSWT],p*p);
53     tr[0] = des->ss[0];
54     tr[1] = M2[0]; /* n int W e^<a,A> */
55   }
56   else
57   { for (i=0; i<des->n; i++)
58     { stdlinks(link,lf,des->ind[i],des->th[i],robscale);
59       ww = SQR(des->w[i])*link[ZDDLL];
60       tr[0] += des->w[i];
61       tr[1] += SQR(des->w[i]);
62       addouter(M2,d_xi(des,i),d_xi(des,i),p,ww);
63     }
64   }
65 
66   memcpy(M12,M2,p*p*sizeof(double));
67   for (i=0; i<p; i++)
68     jacob_solve(&des->xtwx,&M12[i*p]);
69 }
70 
71 /* Compute influence function and estimated derivatives.
72  * Results stored in des->f1.
73  * This assumes weight function is scaled so that W(0)=1.0.
74  */
comp_infl(lf,des)75 double comp_infl(lf,des)
76 lfit *lf;
77 design *des;
78 { unitvec(des->f1,0,des->p);
79   jacob_solve(&des->xtwx,des->f1);
80   return(des->f1[0]);
81 }
82 
comp_vari(lf,des,tr,t0)83 void comp_vari(lf,des,tr,t0)
84 lfit *lf;
85 design *des;
86 double *tr, *t0;
87 { int i, j, k, p;
88   double *M12, *M2;
89   M12 = des->V; M2 = des->P; p = des->p;
90   vmat(lf,des,M12,M2,tr); /* M2 = X^T W^2 V X  tr0=sum(W) tr1=sum(W*W) */
91   tr[2] = m_trace(M12,p);   /* tr (XTWVX)^{-1}(XTW^2VX) */
92 
93   comp_infl(lf,des);
94   for (i=0; i<=lf->mi[MDIM]; i++) t0[i] = des->f1[i];
95 
96 /*
97  * Covariance matrix is M1^{-1} * M2 * M1^{-1}
98  * We compute this using the cholesky decomposition of
99  * M2; premultiplying by M1^{-1} and squaring. This
100  * is more stable than direct computation in near-singular cases.
101  */
102   chol_dec(M2,p);
103   for (i=0; i<p; i++) jacob_solve(&des->xtwx,&M2[i*p]);
104   for (i=0; i<p; i++)
105   { for (j=0; j<p; j++)
106     { M12[i*p+j] = 0;
107       for (k=0; k<p; k++)
108         M12[i*p+j] += M2[k*p+i]*M2[k*p+j]; /* ith column of covariance */
109     }
110   }
111   if ((lf->mi[MTG]==TDEN) && (lf->mi[MLINK]==LIDENT))
112     multmatscal(M12,1/SQR(lf->dp[DSWT]),p*p);
113 }
114 
115 /* local_df computes:
116  *   tr[0] = trace(W)
117  *   tr[1] = trace(W*W)
118  *   tr[2] = trace( M1^{-1} M2 )
119  *   tr[3] = trace( M1^{-1} M3 )
120  *   tr[4] = trace( (M1^{-1} M2)^2 )
121  *   tr[5] = var(theta-hat).
122  */
local_df(lf,des,tr)123 void local_df(lf,des,tr)
124 lfit *lf;
125 design *des;
126 double *tr;
127 { int i, j, p;
128   double *m2, *V, ww, link[LLEN];
129 
130   tr[0] = tr[1] = tr[2] = tr[3] = tr[4] = tr[5] = 0.0;
131   m2 = des->V; V = des->P; p = des->p;
132 
133   vmat(lf,des,m2,V,tr);  /* M = X^T W^2 V X  tr0=sum(W) tr1=sum(W*W) */
134   tr[2] = m_trace(m2,p);   /* tr (XTWVX)^{-1}(XTW^2VX) */
135 
136   unitvec(des->f1,0,p);
137   jacob_solve(&des->xtwx,des->f1);
138   for (i=0; i<p; i++)
139     for (j=0; j<p; j++)
140     { tr[4] += m2[i*p+j]*m2[j*p+i];  /* tr(M^2) */
141       tr[5] += des->f1[i]*V[i*p+j]*des->f1[j]; /* var(thetahat) */
142   }
143   tr[5] = sqrt(tr[5]);
144 
145   setzero(m2,p*p);
146   for (i=0; i<des->n; i++)
147   { stdlinks(link,lf,des->ind[i],des->th[i],robscale);
148     ww = SQR(des->w[i])*des->w[i]*link[ZDDLL];
149     addouter(m2,d_xi(des,i),d_xi(des,i),p,ww);
150   }
151   for (i=0; i<p; i++)
152   { jacob_solve(&des->xtwx,&m2[i*p]);
153     tr[3] += m2[i*(p+1)];
154   }
155 
156   return;
157 }
158