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