Mercurial > repos > vipints > rdiff
diff rDiff/src/locfit/Source/liblocf.c @ 0:0f80a5141704
version 0.3 uploaded
| author | vipints |
|---|---|
| date | Thu, 14 Feb 2013 23:38:36 -0500 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rDiff/src/locfit/Source/liblocf.c Thu Feb 14 23:38:36 2013 -0500 @@ -0,0 +1,4940 @@ +/* + * Copyright 1996-2006 Catherine Loader. + */ + +#include "mex.h" +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + * Integration for hazard rate estimation. The functions in this + * file are used to evaluate + * sum int_0^{Ti} W_i(t,x) A()A()' exp( P() ) dt + * for hazard rate models. + * + * These routines assume the weight function is supported on [-1,1]. + * hasint_sph multiplies by exp(base(lf,i)), which allows estimating + * the baseline in a proportional hazards model, when the covariate + * effect base(lf,i) is known. + * + * TODO: + * hazint_sph, should be able to reduce mint in some cases with + * small integration range. onedint could be used for beta-family + * (RECT,EPAN,BISQ,TRWT) kernels. + * hazint_prod, restrict terms from the sum based on x values. + * I should count obs >= max, and only do that integration once. + */ + +#include "locf.h" + +static double ilim[2*MXDIM], *ff, tmax; +static lfdata *haz_lfd; +static smpar *haz_sp; + +/* + * hrao returns 0 if integration region is empty. + * 1 otherwise. + */ +int haz_sph_int(dfx,cf,h,r1) +double *dfx, *cf, h, *r1; +{ double s, t0, t1, wt, th; + int j, dim, p; + s = 0; p = npar(haz_sp); + dim = haz_lfd->d; + for (j=1; j<dim; j++) s += SQR(dfx[j]/(h*haz_lfd->sca[j])); + if (s>1) return(0); + + setzero(r1,p*p); + t1 = sqrt(1-s)*h*haz_lfd->sca[0]; + t0 = -t1; + if (t0<ilim[0]) t0 = ilim[0]; + if (t1>ilim[dim]) t1 = ilim[dim]; + if (t1>dfx[0]) t1 = dfx[0]; + if (t1<t0) return(0); + +/* Numerical integration by Simpson's rule. + */ + for (j=0; j<=de_mint; j++) + { dfx[0] = t0+(t1-t0)*j/de_mint; + wt = weight(haz_lfd, haz_sp, dfx, NULL, h, 0, 0.0); + fitfun(haz_lfd, haz_sp, dfx,NULL,ff,NULL); + th = innerprod(cf,ff,p); + if (link(haz_sp)==LLOG) th = exp(th); + wt *= 2+2*(j&1)-(j==0)-(j==de_mint); + addouter(r1,ff,ff,p,wt*th); + } + multmatscal(r1,(t1-t0)/(3*de_mint),p*p); + + return(1); +} + +int hazint_sph(t,resp,r1,cf,h) +double *t, *resp, *r1, *cf, h; +{ int i, j, n, p, st; + double dfx[MXDIM], eb, sb; + p = npar(haz_sp); + setzero(resp,p*p); + sb = 0.0; + + n = haz_lfd->n; + for (i=0; i<=n; i++) + { + if (i==n) + { dfx[0] = tmax-t[0]; + for (j=1; j<haz_lfd->d; j++) dfx[j] = 0.0; + eb = exp(sb/n); + } + else + { eb = exp(base(haz_lfd,i)); sb += base(haz_lfd,i); + for (j=0; j<haz_lfd->d; j++) dfx[j] = datum(haz_lfd,j,i)-t[j]; + } + + st = haz_sph_int(dfx,cf,h,r1); + if (st) + for (j=0; j<p*p; j++) resp[j] += eb*r1[j]; + } + return(LF_OK); +} + +int hazint_prod(t,resp,x,cf,h) +double *t, *resp, *x, *cf, h; +{ int d, p, i, j, k, st; + double dfx[MXDIM], t_prev, + hj, hs, ncf[MXDEG], ef, il1; + double prod_wk[MXDIM][2*MXDEG+1], eb, sb; + + p = npar(haz_sp); + d = haz_lfd->d; + setzero(resp,p*p); + hj = hs = h*haz_lfd->sca[0]; + + ncf[0] = cf[0]; + for (i=1; i<=deg(haz_sp); i++) + { ncf[i] = hj*cf[(i-1)*d+1]; hj *= hs; + } + +/* for i=0..n.... + * First we compute prod_wk[j], j=0..d. + * For j=0, this is int_0^T_i (u-t)^k W((u-t)/h) exp(b0*(u-t)) du + * For remaining j, (x(i,j)-x(j))^k Wj exp(bj*(x..-x.)) + * + * Second, we add to the integration (exp(a) incl. in integral) + * with the right factorial denominators. + */ + t_prev = ilim[0]; sb = 0.0; + for (i=0; i<=haz_lfd->n; i++) + { if (i==haz_lfd->n) + { dfx[0] = tmax-t[0]; + for (j=1; j<d; j++) dfx[j] = 0.0; + eb = exp(sb/haz_lfd->n); + } + else + { eb = exp(base(haz_lfd,i)); sb += base(haz_lfd,i); + for (j=0; j<d; j++) dfx[j] = datum(haz_lfd,j,i)-t[j]; + } + + if (dfx[0]>ilim[0]) /* else it doesn't contribute */ + { +/* time integral */ + il1 = (dfx[0]>ilim[d]) ? ilim[d] : dfx[0]; + if (il1 != t_prev) /* don't repeat! */ + { st = onedint(haz_sp,ncf,ilim[0]/hs,il1/hs,prod_wk[0]); + if (st>0) return(st); + hj = eb; + for (j=0; j<=2*deg(haz_sp); j++) + { hj *= hs; + prod_wk[0][j] *= hj; + } + t_prev = il1; + } + +/* covariate terms */ + for (j=1; j<d; j++) + { + ef = 0.0; + for (k=deg(haz_sp); k>0; k--) ef = (ef+dfx[j])*cf[1+(k-1)*d+j]; + ef = exp(ef); + prod_wk[j][0] = ef * W(dfx[j]/(h*haz_lfd->sca[j]),ker(haz_sp)); + for (k=1; k<=2*deg(haz_sp); k++) + prod_wk[j][k] = prod_wk[j][k-1] * dfx[j]; + } + +/* add to the integration. */ + prodintresp(resp,prod_wk,d,deg(haz_sp),p); + } /* if dfx0 > ilim0 */ + } /* n loop */ + +/* symmetrize */ + for (k=0; k<p; k++) + for (j=k; j<p; j++) + resp[j*p+k] = resp[k*p+j]; + return(LF_OK); +} + +int hazint(t,resp,resp1,cf,h) +double *t, *resp, *resp1, *cf, h; +{ if (haz_lfd->d==1) return(hazint_prod(t,resp,resp1,cf,h)); + if (kt(haz_sp)==KPROD) return(hazint_prod(t,resp,resp1,cf,h)); + + return(hazint_sph(t,resp,resp1,cf,h)); +} + +void haz_init(lfd,des,sp,il) +lfdata *lfd; +design *des; +smpar *sp; +double *il; +{ int i; + + haz_lfd = lfd; + haz_sp = sp; + + tmax = datum(lfd,0,0); + for (i=1; i<lfd->n; i++) tmax = MAX(tmax,datum(lfd,0,i)); + ff = des->xtwx.wk; + for (i=0; i<2*lfd->d; i++) ilim[i] = il[i]; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + * + * Routines for one-dimensional numerical integration + * in density estimation. The entry point is + * + * onedint(cf,mi,l0,l1,resp) + * + * which evaluates int W(u)u^j exp( P(u) ), j=0..2*deg. + * P(u) = cf[0] + cf[1]u + cf[2]u^2/2 + ... + cf[deg]u^deg/deg! + * l0 and l1 are the integration limits. + * The results are returned through the vector resp. + * + */ + +#include "locf.h" + +static int debug; + +int exbctay(b,c,n,z) /* n-term taylor series of e^(bx+cx^2) */ +double b, c, *z; +int n; +{ double ec[20]; + int i, j; + z[0] = 1; + for (i=1; i<=n; i++) z[i] = z[i-1]*b/i; + if (c==0.0) return(n); + if (n>=40) + { WARN(("exbctay limit to n<40")); + n = 39; + } + ec[0] = 1; + for (i=1; 2*i<=n; i++) ec[i] = ec[i-1]*c/i; + for (i=n; i>1; i--) + for (j=1; 2*j<=i; j++) + z[i] += ec[j]*z[i-2*j]; + return(n); +} + +double explinjtay(l0,l1,j,cf) +/* int_l0^l1 x^j e^(a+bx+cx^2); exbctay aroud l1 */ +double l0, l1, *cf; +int j; +{ double tc[40], f, s; + int k, n; + if ((l0!=0.0) | (l1!=1.0)) WARN(("explinjtay: invalid l0, l1")); + n = exbctay(cf[1]+2*cf[2]*l1,cf[2],20,tc); + s = tc[0]/(j+1); + f = 1/(j+1); + for (k=1; k<=n; k++) + { f *= -k/(j+k+1.0); + s += tc[k]*f; + } + return(f); +} + +void explint1(l0,l1,cf,I,p) /* int x^j exp(a+bx); j=0..p-1 */ +double l0, l1, *cf, *I; +int p; +{ double y0, y1, f; + int j, k, k1; + y0 = mut_exp(cf[0]+l0*cf[1]); + y1 = mut_exp(cf[0]+l1*cf[1]); + if (p<2*fabs(cf[1])) k = p; else k = (int)fabs(cf[1]); + + if (k>0) + { I[0] = (y1-y0)/cf[1]; + for (j=1; j<k; j++) /* forward steps for small j */ + { y1 *= l1; y0 *= l0; + I[j] = (y1-y0-j*I[j-1])/cf[1]; + } + if (k==p) return; + y1 *= l1; y0 *= l0; + } + + f = 1; k1 = k; + while ((k<50) && (f>1.0e-8)) /* initially Ik = diff(x^{k+1}e^{a+bx}) */ + { y1 *= l1; y0 *= l0; + I[k] = y1-y0; + if (k>=p) f *= fabs(cf[1])/(k+1); + k++; + } + if (k==50) WARN(("explint1: want k>50")); + I[k] = 0.0; + for (j=k-1; j>=k1; j--) /* now do back step recursion */ + I[j] = (I[j]-cf[1]*I[j+1])/(j+1); +} + +void explintyl(l0,l1,cf,I,p) /* small c, use taylor series and explint1 */ +double l0, l1, *cf, *I; +int p; +{ int i; + double c; + explint1(l0,l1,cf,I,p+8); + c = cf[2]; + for (i=0; i<p; i++) + I[i] = (((I[i+8]*c/4+I[i+6])*c/3+I[i+4])*c/2+I[i+2])*c+I[i]; +} + +void solvetrid(X,y,m) +double *X, *y; +int m; +{ int i; + double s; + for (i=1; i<m; i++) + { s = X[3*i]/X[3*i-2]; + X[3*i] = 0; X[3*i+1] -= s*X[3*i-1]; + y[i] -= s*y[i-1]; + } + for (i=m-2; i>=0; i--) + { s = X[3*i+2]/X[3*i+4]; + X[3*i+2] = 0; + y[i] -= s*y[i+1]; + } + for (i=0; i<m; i++) y[i] /= X[3*i+1]; +} + +void initi0i1(I,cf,y0,y1,l0,l1) +double *I, *cf, y0, y1, l0, l1; +{ double a0, a1, c, d, bi; + d = -cf[1]/(2*cf[2]); c = sqrt(2*fabs(cf[2])); + a0 = c*(l0-d); a1 = c*(l1-d); + if (cf[2]<0) + { bi = mut_exp(cf[0]+cf[1]*d+cf[2]*d*d)/c; + if (a0>0) + { if (a0>6) I[0] = (y0*ptail(-a0)-y1*ptail(-a1))/c; + else I[0] = S2PI*(mut_pnorm(-a0)-mut_pnorm(-a1))*bi; + } + else + { if (a1< -6) I[0] = (y1*ptail(a1)-y0*ptail(a0))/c; + else I[0] = S2PI*(mut_pnorm(a1)-mut_pnorm(a0))*bi; + } + } + else + I[0] = (y1*mut_daws(a1)-y0*mut_daws(a0))/c; + I[1] = (y1-y0)/(2*cf[2])+d*I[0]; +} + +void explinsid(l0,l1,cf,I,p) /* large b; don't use fwd recursion */ +double l0, l1, *cf, *I; +int p; +{ int k, k0, k1, k2; + double y0, y1, Z[150]; +if (debug) mut_printf("side: %8.5f %8.5f %8.5f limt %8.5f %8.5f p %2d\n",cf[0],cf[1],cf[2],l0,l1,p); + + k0 = 2; + k1 = (int)(fabs(cf[1])+fabs(2*cf[2])); + if (k1<2) k1 = 2; + if (k1>p+20) k1 = p+20; + k2 = p+20; + +if (k2>50) { mut_printf("onedint: k2 warning\n"); k2 = 50; } + if (debug) mut_printf("k0 %2d k1 %2d k2 %2d p %2d\n",k0,k1,k2,p); + + y0 = mut_exp(cf[0]+l0*(cf[1]+l0*cf[2])); + y1 = mut_exp(cf[0]+l1*(cf[1]+l1*cf[2])); + initi0i1(I,cf,y0,y1,l0,l1); +if (debug) mut_printf("i0 %8.5f i1 %8.5f\n",I[0],I[1]); + + y1 *= l1; y0 *= l0; /* should be x^(k1)*exp(..) */ + if (k0<k1) /* center steps; initially x^k*exp(...) */ + for (k=k0; k<k1; k++) + { y1 *= l1; y0 *= l0; + I[k] = y1-y0; + Z[3*k] = k; Z[3*k+1] = cf[1]; Z[3*k+2] = 2*cf[2]; + } + + y1 *= l1; y0 *= l0; /* should be x^(k1)*exp(..) */ +if (debug) mut_printf("k1 %2d y0 %8.5f y1 %8.5f\n",k1,y0,y1); + for (k=k1; k<k2; k++) + { y1 *= l1; y0 *= l0; + I[k] = y1-y0; + } + I[k2] = I[k2+1] = 0.0; + for (k=k2-1; k>=k1; k--) + I[k] = (I[k]-cf[1]*I[k+1]-2*cf[2]*I[k+2])/(k+1); + + if (k0<k1) + { I[k0] -= k0*I[k0-1]; + I[k1-1] -= 2*cf[2]*I[k1]; + Z[3*k0] = Z[3*k1-1] = 0; + solvetrid(&Z[3*k0],&I[k0],k1-k0); + } +if (debug) +{ mut_printf("explinsid:\n"); + for (k=0; k<p; k++) mut_printf(" %8.5f\n",I[k]); +} +} + +void explinbkr(l0,l1,cf,I,p) /* small b,c; use back recursion */ +double l0, l1, *cf, *I; +int p; +{ int k, km; + double y0, y1; + y0 = mut_exp(cf[0]+l0*(cf[1]+cf[2]*l0)); + y1 = mut_exp(cf[0]+l1*(cf[1]+cf[2]*l1)); + km = p+10; + for (k=0; k<=km; k++) + { y1 *= l1; y0 *= l0; + I[k] = y1-y0; + } + I[km+1] = I[km+2] = 0; + for (k=km; k>=0; k--) + I[k] = (I[k]-cf[1]*I[k+1]-2*cf[2]*I[k+2])/(k+1); +} + +void explinfbk0(l0,l1,cf,I,p) /* fwd and bac recur; b=0; c<0 */ +double l0, l1, *cf, *I; +int p; +{ double y0, y1, f1, f2, f, ml2; + int k, ks; + + y0 = mut_exp(cf[0]+l0*l0*cf[2]); + y1 = mut_exp(cf[0]+l1*l1*cf[2]); + initi0i1(I,cf,y0,y1,l0,l1); + + ml2 = MAX(l0*l0,l1*l1); + ks = 1+(int)(2*fabs(cf[2])*ml2); + if (ks<2) ks = 2; + if (ks>p-3) ks = p; + + /* forward recursion for k < ks */ + for (k=2; k<ks; k++) + { y1 *= l1; y0 *= l0; + I[k] = (y1-y0-(k-1)*I[k-2])/(2*cf[2]); + } + if (ks==p) return; + + y1 *= l1*l1; y0 *= l0*l0; + for (k=ks; k<p; k++) /* set I[k] = x^{k+1}e^(a+cx^2) | {l0,l1} */ + { y1 *= l1; y0 *= l0; + I[k] = y1-y0; + } + + /* initialize I[p-2] and I[p-1] */ + f1 = 1.0/p; f2 = 1.0/(p-1); + I[p-1] *= f1; I[p-2] *= f2; + k = p; f = 1.0; + while (f>1.0e-8) + { y1 *= l1; y0 *= l0; + if ((k-p)%2==0) /* add to I[p-2] */ + { f2 *= -2*cf[2]/(k+1); + I[p-2] += (y1-y0)*f2; + } + else /* add to I[p-1] */ + { f1 *= -2*cf[2]/(k+1); + I[p-1] += (y1-y0)*f1; + f *= 2*fabs(cf[2])*ml2/(k+1); + } + k++; + } + + /* use back recursion for I[ks..(p-3)] */ + for (k=p-3; k>=ks; k--) + I[k] = (I[k]-2*cf[2]*I[k+2])/(k+1); +} + +void explinfbk(l0,l1,cf,I,p) /* fwd and bac recur; b not too large */ +double l0, l1, *cf, *I; +int p; +{ double y0, y1; + int k, ks, km; + + y0 = mut_exp(cf[0]+l0*(cf[1]+l0*cf[2])); + y1 = mut_exp(cf[0]+l1*(cf[1]+l1*cf[2])); + initi0i1(I,cf,y0,y1,l0,l1); + + ks = (int)(3*fabs(cf[2])); + if (ks<3) ks = 3; + if (ks>0.75*p) ks = p; /* stretch the forward recurs as far as poss. */ + /* forward recursion for k < ks */ + for (k=2; k<ks; k++) + { y1 *= l1; y0 *= l0; + I[k] = (y1-y0-cf[1]*I[k-1]-(k-1)*I[k-2])/(2*cf[2]); + } + if (ks==p) return; + + km = p+15; + y1 *= l1*l1; y0 *= l0*l0; + for (k=ks; k<=km; k++) + { y1 *= l1; y0 *= l0; + I[k] = y1-y0; + } + I[km+1] = I[km+2] = 0.0; + for (k=km; k>=ks; k--) + I[k] = (I[k]-cf[1]*I[k+1]-2*cf[2]*I[k+2])/(k+1); +} + +void recent(I,resp,wt,p,s,x) +double *I, *resp, *wt, x; +int p, s; +{ int i, j; + + /* first, use W taylor series I -> resp */ + for (i=0; i<=p; i++) + { resp[i] = 0.0; + for (j=0; j<s; j++) resp[i] += wt[j]*I[i+j]; + } + + /* now, recenter x -> 0 */ + if (x==0) return; + for (j=0; j<=p; j++) for (i=p; i>j; i--) resp[i] += x*resp[i-1]; +} + +void recurint(l0,l2,cf,resp,p,ker) +double l0, l2, *cf, *resp; +int p, ker; +{ int i, s; + double l1, d0, d1, d2, dl, z0, z1, z2, wt[20], ncf[3], I[50], r1[5], r2[5]; +if (debug) mut_printf("\nrecurint: %8.5f %8.5f %8.5f %8.5f %8.5f\n",cf[0],cf[1],cf[2],l0,l2); + + if (cf[2]==0) /* go straight to explint1 */ + { s = wtaylor(wt,0.0,ker); +if (debug) mut_printf("case 1\n"); + explint1(l0,l2,cf,I,p+s); + recent(I,resp,wt,p,s,0.0); + return; + } + + dl = l2-l0; + d0 = cf[1]+2*l0*cf[2]; + d2 = cf[1]+2*l2*cf[2]; + z0 = cf[0]+l0*(cf[1]+l0*cf[2]); + z2 = cf[0]+l2*(cf[1]+l2*cf[2]); + + if ((fabs(cf[1]*dl)<1) && (fabs(cf[2]*dl*dl)<1)) + { ncf[0] = z0; ncf[1] = d0; ncf[2] = cf[2]; +if (debug) mut_printf("case 2\n"); + s = wtaylor(wt,l0,ker); + explinbkr(0.0,dl,ncf,I,p+s); + recent(I,resp,wt,p,s,l0); + return; + } + + if (fabs(cf[2]*dl*dl)<0.001) /* small c, use explint1+tay.ser */ + { ncf[0] = z0; ncf[1] = d0; ncf[2] = cf[2]; +if (debug) mut_printf("case small c\n"); + s = wtaylor(wt,l0,ker); + explintyl(0.0,l2-l0,ncf,I,p+s); + recent(I,resp,wt,p,s,l0); + return; + } + + if (d0*d2<=0) /* max/min in [l0,l2] */ + { l1 = -cf[1]/(2*cf[2]); + z1 = cf[0]+l1*(cf[1]+l1*cf[2]); + d1 = 0.0; + if (cf[2]<0) /* peak, integrate around l1 */ + { s = wtaylor(wt,l1,ker); + ncf[0] = z1; ncf[1] = 0.0; ncf[2] = cf[2]; +if (debug) mut_printf("case peak p %2d s %2d\n",p,s); + explinfbk0(l0-l1,l2-l1,ncf,I,p+s); + recent(I,resp,wt,p,s,l1); + return; + } + } + + if ((d0-2*cf[2]*dl)*(d2+2*cf[2]*dl)<0) /* max/min is close to [l0,l2] */ + { l1 = -cf[1]/(2*cf[2]); + z1 = cf[0]+l1*(cf[1]+l1*cf[2]); + if (l1<l0) { l1 = l0; z1 = z0; } + if (l1>l2) { l1 = l2; z1 = z2; } + + if ((z1>=z0) & (z1>=z2)) /* peak; integrate around l1 */ + { s = wtaylor(wt,l1,ker); +if (debug) mut_printf("case 4\n"); + d1 = cf[1]+2*l1*cf[2]; + ncf[0] = z1; ncf[1] = d1; ncf[2] = cf[2]; + explinfbk(l0-l1,l2-l1,ncf,I,p+s); + recent(I,resp,wt,p,s,l1); + return; + } + + /* trough; integrate [l0,l1] and [l1,l2] */ + for (i=0; i<=p; i++) r1[i] = r2[i] = 0.0; + if (l0<l1) + { s = wtaylor(wt,l0,ker); +if (debug) mut_printf("case 5\n"); + ncf[0] = z0; ncf[1] = d0; ncf[2] = cf[2]; + explinfbk(0.0,l1-l0,ncf,I,p+s); + recent(I,r1,wt,p,s,l0); + } + if (l1<l2) + { s = wtaylor(wt,l2,ker); +if (debug) mut_printf("case 6\n"); + ncf[0] = z2; ncf[1] = d2; ncf[2] = cf[2]; + explinfbk(l1-l2,0.0,ncf,I,p+s); + recent(I,r2,wt,p,s,l2); + } + for (i=0; i<=p; i++) resp[i] = r1[i]+r2[i]; + return; + } + + /* Now, quadratic is monotone on [l0,l2]; big b; moderate c */ + if (z2>z0+3) /* steep increase, expand around l2 */ + { s = wtaylor(wt,l2,ker); +if (debug) mut_printf("case 7\n"); + + + ncf[0] = z2; ncf[1] = d2; ncf[2] = cf[2]; + explinsid(l0-l2,0.0,ncf,I,p+s); + recent(I,resp,wt,p,s,l2); +if (debug) mut_printf("7 resp: %8.5f %8.5f %8.5f %8.5f\n",resp[0],resp[1],resp[2],resp[3]); + return; + } + + /* bias towards expansion around l0, because it's often 0 */ +if (debug) mut_printf("case 8\n"); + s = wtaylor(wt,l0,ker); + ncf[0] = z0; ncf[1] = d0; ncf[2] = cf[2]; + explinsid(0.0,l2-l0,ncf,I,p+s); + recent(I,resp,wt,p,s,l0); + return; +} + +int onedexpl(cf,deg,resp) +double *cf, *resp; +int deg; +{ int i; + double f0, fr, fl; + if (deg>=2) LERR(("onedexpl only valid for deg=0,1")); + if (fabs(cf[1])>=EFACT) return(LF_BADP); + + f0 = exp(cf[0]); fl = fr = 1.0; + for (i=0; i<=2*deg; i++) + { f0 *= i+1; + fl /=-(EFACT+cf[1]); + fr /= EFACT-cf[1]; + resp[i] = f0*(fr-fl); + } + return(LF_OK); +} + +int onedgaus(cf,deg,resp) +double *cf, *resp; +int deg; +{ int i; + double f0, mu, s2; + if (deg==3) + { LERR(("onedgaus only valid for deg=0,1,2")); + return(LF_ERR); + } + if (2*cf[2]>=GFACT*GFACT) return(LF_BADP); + + s2 = 1/(GFACT*GFACT-2*cf[2]); + mu = cf[1]*s2; + resp[0] = 1.0; + if (deg>=1) + { resp[1] = mu; + resp[2] = s2+mu*mu; + if (deg==2) + { resp[3] = mu*(3*s2+mu*mu); + resp[4] = 3*s2*s2 + mu*mu*(6*s2+mu*mu); + } + } + f0 = S2PI * exp(cf[0]+mu*mu/(2*s2))*sqrt(s2); + for (i=0; i<=2*deg; i++) resp[i] *= f0; + return(LF_OK); +} + +int onedint(sp,cf,l0,l1,resp) /* int W(u)u^j exp(..), j=0..2*deg */ +smpar *sp; +double *cf, l0, l1, *resp; +{ double u, uj, y, ncf[4], rr[5]; + int i, j; + +if (debug) mut_printf("onedint: %f %f %f %f %f\n",cf[0],cf[1],cf[2],l0,l1); + + if (deg(sp)<=2) + { for (i=0; i<3; i++) ncf[i] = (i>deg(sp)) ? 0.0 : cf[i]; + ncf[2] /= 2; + + if (ker(sp)==WEXPL) return(onedexpl(ncf,deg(sp),resp)); + if (ker(sp)==WGAUS) return(onedgaus(ncf,deg(sp),resp)); + + if (l1>0) + recurint(MAX(l0,0.0),l1,ncf,resp,2*deg(sp),ker(sp)); + else for (i=0; i<=2*deg(sp); i++) resp[i] = 0; + + if (l0<0) + { ncf[1] = -ncf[1]; + l0 = -l0; l1 = -l1; + recurint(MAX(l1,0.0),l0,ncf,rr,2*deg(sp),ker(sp)); + } + else for (i=0; i<=2*deg(sp); i++) rr[i] = 0.0; + + for (i=0; i<=2*deg(sp); i++) + resp[i] += (i%2==0) ? rr[i] : -rr[i]; + + return(LF_OK); + } + + /* For degree >= 3, we use Simpson's rule. */ + for (j=0; j<=2*deg(sp); j++) resp[j] = 0.0; + for (i=0; i<=de_mint; i++) + { u = l0+(l1-l0)*i/de_mint; + y = cf[0]; uj = 1; + for (j=1; j<=deg(sp); j++) + { uj *= u; + y += cf[j]*uj/fact[j]; + } + y = (4-2*(i%2==0)-(i==0)-(i==de_mint)) * + W(fabs(u),ker(sp))*exp(MIN(y,300.0)); + for (j=0; j<=2*deg(sp); j++) + { resp[j] += y; + y *= u; + } + } + for (j=0; j<=2*deg(sp); j++) resp[j] = resp[j]*(l1-l0)/(3*de_mint); + return(LF_OK); +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +extern int lf_status; +static double u[MXDIM], ilim[2*MXDIM], *ff, hh, *cff; +static lfdata *den_lfd; +static design *den_des; +static smpar *den_sp; +int fact[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800}; +int de_mint = 20; +int de_itype = IDEFA; +int de_renorm= 0; + +int multint(), prodint(), gausint(), mlinint(); + +#define NITYPE 7 +static char *itype[NITYPE] = { "default", "multi", "product", "mlinear", + "hazard", "sphere", "monte" }; +static int ivals[NITYPE] = + { IDEFA, IMULT, IPROD, IMLIN, IHAZD, ISPHR, IMONT }; +int deitype(char *z) +{ return(pmatch(z, itype, ivals, NITYPE, IDEFA)); +} + +void prresp(coef,resp,p) +double *coef, *resp; +int p; +{ int i, j; + mut_printf("Coefficients:\n"); + for (i=0; i<p; i++) mut_printf("%8.5f ",coef[i]); + mut_printf("\n"); + mut_printf("Response matrix:\n"); + for (i=0; i<p; i++) + { for (j=0; j<p; j++) mut_printf("%9.6f, ",resp[i+j*p]); + mut_printf("\n"); + } +} + +int mif(u,d,resp,M) +double *u, *resp, *M; +int d; +{ double wt; + int i, j, p; + + p = den_des->p; + wt = weight(den_lfd, den_sp, u, NULL, hh, 0, 0.0); + if (wt==0) + { setzero(resp,p*p); + return(p*p); + } + + fitfun(den_lfd, den_sp, u,NULL,ff,NULL); + if (link(den_sp)==LLOG) + wt *= mut_exp(innerprod(ff,cff,p)); + for (i=0; i<p; i++) + for (j=0; j<p; j++) + resp[i*p+j] = wt*ff[i]*ff[j]; + return(p*p); +} + +int multint(t,resp1,resp2,cf,h) +double *t, *resp1, *resp2, *cf, h; +{ int d, i, mg[MXDIM]; + + if (ker(den_sp)==WGAUS) return(gausint(t,resp1,resp2,cf,h,den_lfd->sca)); + + d = den_lfd->d; + for (i=0; i<d; i++) mg[i] = de_mint; + + hh = h; + cff= cf; + simpsonm(mif,ilim,&ilim[d],d,resp1,mg,resp2); + return(LF_OK); +} + +int mlinint(t,resp1,resp2,cf,h) +double *t, *resp1, *resp2, *cf, h; +{ + double hd, nb, wt, wu, g[4], w0, w1, v, *sca; + int d, p, i, j, jmax, k, l, z, jj[2]; + + d = den_lfd->d; p = den_des->p; sca = den_lfd->sca; + hd = 1; + for (i=0; i<d; i++) hd *= h*sca[i]; + + if (link(den_sp)==LIDENT) + { setzero(resp1,p*p); + resp1[0] = wint(d,NULL,0,ker(den_sp))*hd; + if (deg(den_sp)==0) return(LF_OK); + jj[0] = 2; w0 = wint(d,jj,1,ker(den_sp))*hd*h*h; + for (i=0; i<d; i++) resp1[(i+1)*p+i+1] = w0*sca[i]*sca[i]; + if (deg(den_sp)==1) return(LF_OK); + for (i=0; i<d; i++) + { j = p-(d-i)*(d-i+1)/2; + resp1[j] = resp1[p*j] = w0*sca[i]*sca[i]/2; + } + if (d>1) + { jj[1] = 2; + w0 = wint(d,jj,2,ker(den_sp)) * hd*h*h*h*h; + } + jj[0] = 4; + w1 = wint(d,jj,1,ker(den_sp)) * hd*h*h*h*h/4; + z = d+1; + for (i=0; i<d; i++) + { k = p-(d-i)*(d-i+1)/2; + for (j=i; j<d; j++) + { l = p-(d-j)*(d-j+1)/2; + if (i==j) resp1[z*p+z] = w1*SQR(sca[i])*SQR(sca[i]); + else + { resp1[z*p+z] = w0*SQR(sca[i])*SQR(sca[j]); + resp1[k*p+l] = resp1[k+p*l] = w0/4*SQR(sca[i])*SQR(sca[j]); + } + z++; + } } + return(LF_OK); + } + switch(deg(den_sp)) + { case 0: + resp1[0] = mut_exp(cf[0])*wint(d,NULL,0,ker(den_sp))*hd; + return(LF_OK); + case 1: + nb = 0.0; + for (i=1; i<=d; i++) + { v = h*cf[i]*sca[i-1]; + nb += v*v; + } + if (ker(den_sp)==WGAUS) + { w0 = 1/(GFACT*GFACT); + g[0] = mut_exp(cf[0]+w0*nb/2+d*log(S2PI/2.5)); + g[1] = g[3] = g[0]*w0; + g[2] = g[0]*w0*w0; + } + else + { wt = wu = mut_exp(cf[0]); + w0 = wint(d,NULL,0,ker(den_sp)); g[0] = wt*w0; + g[1] = g[2] = g[3] = 0.0; + j = 0; jmax = (d+2)*de_mint; + while ((j<jmax) && (wt*w0/g[0]>1.0e-8)) + { j++; + jj[0] = 2*j; w0 = wint(d,jj,1,ker(den_sp)); + if (d==1) g[3] += wt * w0; + else + { jj[0] = 2; jj[1] = 2*j-2; w1 = wint(d,jj,2,ker(den_sp)); + g[3] += wt*w1; + g[2] += wu*(w0-w1); + } + wt /= (2*j-1.0); g[1] += wt*w0; + wt *= nb/(2*j); g[0] += wt*w0; + wu /= (2*j-1.0)*(2*j); + if (j>1) wu *= nb; + } + if (j==jmax) WARN(("mlinint: series not converged")); + } + g[0] *= hd; g[1] *= hd; + g[2] *= hd; g[3] *= hd; + resp1[0] = g[0]; + for (i=1; i<=d; i++) + { resp1[i] = resp1[(d+1)*i] = cf[i]*SQR(h*sca[i-1])*g[1]; + for (j=1; j<=d; j++) + { resp1[(d+1)*i+j] = (i==j) ? g[3]*SQR(h*sca[i-1]) : 0; + resp1[(d+1)*i+j] += g[2]*SQR(h*h*sca[i-1]*sca[j-1])*cf[i]*cf[j]; + } + } + return(LF_OK); + } + LERR(("mlinint: deg=0,1 only")); + return(LF_ERR); +} + +void prodintresp(resp,prod_wk,dim,deg,p) +double *resp, prod_wk[MXDIM][2*MXDEG+1]; +int dim, deg, p; +{ double prod; + int i, j, k, j1, k1; + + prod = 1.0; + for (i=0; i<dim; i++) prod *= prod_wk[i][0]; + resp[0] += prod; + if (deg==0) return; + + for (j1=1; j1<=deg; j1++) + { for (j=0; j<dim; j++) + { prod = 1.0; + for (i=0; i<dim; i++) prod *= prod_wk[i][j1*(j==i)]; + prod /= fact[j1]; + resp[1 + (j1-1)*dim +j] += prod; + } + } + + for (k1=1; k1<=deg; k1++) + for (j1=k1; j1<=deg; j1++) + { for (k=0; k<dim; k++) + for (j=0; j<dim; j++) + { prod = 1.0; + for (i=0; i<dim; i++) prod *= prod_wk[i][k1*(k==i) + j1*(j==i)]; + prod /= fact[k1]*fact[j1]; + resp[ (1+(k1-1)*dim+k)*p + 1+(j1-1)*dim+j] += prod; + } + } +} + +int prodint(t,resp,resp2,coef,h) +double *t, *resp, *resp2, *coef, h; +{ int dim, p, i, j, k, st; + double cf[MXDEG+1], hj, hs, prod_wk[MXDIM][2*MXDEG+1]; + + dim = den_lfd->d; + p = den_des->p; + for (i=0; i<p*p; i++) resp[i] = 0.0; + cf[0] = coef[0]; + +/* compute the one dimensional terms + */ + for (i=0; i<dim; i++) + { hj = 1; hs = h*den_lfd->sca[i]; + for (j=0; j<deg(den_sp); j++) + { hj *= hs; + cf[j+1] = hj*coef[ j*dim+i+1 ]; + } + st = onedint(den_sp,cf,ilim[i]/hs,ilim[i+dim]/hs,prod_wk[i]); + if (st==LF_BADP) return(st); + hj = 1; + for (j=0; j<=2*deg(den_sp); j++) + { hj *= hs; + prod_wk[i][j] *= hj; + } + cf[0] = 0.0; /* so we only include it once, when d>=2 */ + } + +/* transfer to the resp array + */ + prodintresp(resp,prod_wk,dim,deg(den_sp),p); + +/* Symmetrize. +*/ + for (k=0; k<p; k++) + for (j=k; j<p; j++) + resp[j*p+k] = resp[k*p+j]; + + return(st); +} + +int gausint(t,resp,C,cf,h,sca) +double *t, *resp, *C, *cf, h, *sca; +{ double nb, det, z, *P; + int d, p, i, j, k, l, m1, m2, f; + d = den_lfd->d; p = den_des->p; + m1 = d+1; nb = 0; + P = &C[d*d]; + resp[0] = 1; + for (i=0; i<d; i++) + { C[i*d+i] = SQR(GFACT/(h*sca[i]))-cf[m1++]; + for (j=i+1; j<d; j++) C[i*d+j] = C[j*d+i] = -cf[m1++]; + } + eig_dec(C,P,d); + det = 1; + for (i=1; i<=d; i++) + { det *= C[(i-1)*(d+1)]; + if (det <= 0) return(LF_BADP); + resp[i] = cf[i]; + for (j=1; j<=d; j++) resp[j+i*p] = 0; + resp[i+i*p] = 1; + svdsolve(&resp[i*p+1],u,P,C,P,d,0.0); + } + svdsolve(&resp[1],u,P,C,P,d,0.0); + det = sqrt(det); + for (i=1; i<=d; i++) + { nb += cf[i]*resp[i]; + resp[i*p] = resp[i]; + for (j=1; j<=d; j++) + resp[i+p*j] += resp[i]*resp[j]; + } + m1 = d; + for (i=1; i<=d; i++) + for (j=i; j<=d; j++) + { m1++; f = 1+(i==j); + resp[m1] = resp[m1*p] = resp[i*p+j]/f; + m2 = d; + for (k=1; k<=d; k++) + { resp[m1+k*p] = resp[k+m1*p] = + ( resp[i]*resp[j*p+k] + resp[j]*resp[i*p+k] + + resp[k]*resp[i*p+j] - 2*resp[i]*resp[j]*resp[k] )/f; + for (l=k; l<=d; l++) + { m2++; f = (1+(i==j))*(1+(k==l)); + resp[m1+m2*p] = resp[m2+m1*p] = ( resp[i+j*p]*resp[k+l*p] + + resp[i+k*p]*resp[j+l*p] + resp[i+l*p]*resp[j+k*p] + - 2*resp[i]*resp[j]*resp[k]*resp[l] )/f; + } } } + z = mut_exp(d*0.918938533+cf[0]+nb/2)/det; + multmatscal(resp,z,p*p); + return(LF_OK); +} + +int likeden(coef, lk0, f1, A) +double *coef, *lk0, *f1, *A; +{ double lk, r; + int i, j, p, rstat; + + lf_status = LF_OK; + p = den_des->p; + if ((link(den_sp)==LIDENT) && (coef[0] != 0.0)) return(NR_BREAK); + lf_status = (den_des->itype)(den_des->xev,A,den_des->xtwx.Q,coef,den_des->h); + if (lf_error) lf_status = LF_ERR; + if (lf_status==LF_BADP) + { *lk0 = -1.0e300; + return(NR_REDUCE); + } + if (lf_status!=LF_OK) return(NR_BREAK); + if (lf_debug>2) prresp(coef,A,p); + + den_des->xtwx.p = p; + rstat = NR_OK; + switch(link(den_sp)) + { case LLOG: + r = den_des->ss[0]/A[0]; + coef[0] += log(r); + multmatscal(A,r,p*p); + A[0] = den_des->ss[0]; + lk = -A[0]; + if (fabs(coef[0]) > 700) + { lf_status = LF_OOB; + rstat = NR_REDUCE; + } + for (i=0; i<p; i++) + { lk += coef[i]*den_des->ss[i]; + f1[i] = den_des->ss[i]-A[i]; + } + break; + case LIDENT: + lk = 0.0; + for (i=0; i<p; i++) + { f1[i] = den_des->ss[i]; + for (j=0; j<p; j++) + den_des->res[i] -= A[i*p+j]*coef[j]; + } + break; + } + *lk0 = den_des->llk = lk; + + return(rstat); +} + +int inre(x,bound,d) +double *x, *bound; +int d; +{ int i, z; + z = 1; + for (i=0; i<d; i++) + if (bound[i]<bound[i+d]) + z &= (x[i]>=bound[i]) & (x[i]<=bound[i+d]); + return(z); +} + +int setintlimits(lfd, x, h, ang, lset) +lfdata *lfd; +int *ang, *lset; +double *x, h; +{ int d, i; + d = lfd->d; + *ang = *lset = 0; + for (i=0; i<d; i++) + { if (lfd->sty[i]==STANGL) + { ilim[i+d] = ((h<2) ? 2*asin(h/2) : PI)*lfd->sca[i]; + ilim[i] = -ilim[i+d]; + *ang = 1; + } + else + { ilim[i+d] = h*lfd->sca[i]; + ilim[i] = -ilim[i+d]; + + if (lfd->sty[i]==STLEFT) { ilim[i+d] = 0; *lset = 1; } + if (lfd->sty[i]==STRIGH) { ilim[i] = 0; *lset = 1; } + + if (lfd->xl[i]<lfd->xl[i+d]) /* user limits for this variable */ + { if (lfd->xl[i]-x[i]> ilim[i]) + { ilim[i] = lfd->xl[i]-x[i]; *lset=1; } + if (lfd->xl[i+d]-x[i]< ilim[i+d]) + { ilim[i+d] = lfd->xl[i+d]-x[i]; *lset=1; } + } + } + if (ilim[i]==ilim[i+d]) return(LF_DEMP); /* empty integration */ + } + return(LF_OK); +} + +int selectintmeth(itype,lset,ang) +int itype, lset, ang; +{ + if (itype==IDEFA) /* select the default method */ + { if (fam(den_sp)==THAZ) + { if (ang) return(IDEFA); + return( IHAZD ); + } + + if (ubas(den_sp)) return(IMULT); + + if (ang) return(IMULT); + + if (iscompact(ker(den_sp))) + { if (kt(den_sp)==KPROD) return(IPROD); + if (lset) + return( (den_lfd->d==1) ? IPROD : IMULT ); + if (deg(den_sp)<=1) return(IMLIN); + if (den_lfd->d==1) return(IPROD); + return(IMULT); + } + + if (ker(den_sp)==WGAUS) + { if (lset) WARN(("Integration for Gaussian weights ignores limits")); + if ((den_lfd->d==1)|(kt(den_sp)==KPROD)) return(IPROD); + if (deg(den_sp)<=1) return(IMLIN); + if (deg(den_sp)==2) return(IMULT); + } + + return(IDEFA); + } + + /* user provided an integration method, check it is valid */ + + if (fam(den_sp)==THAZ) + { if (ang) return(INVLD); + if (!iscompact(ker(den_sp))) return(INVLD); + return( ((kt(den_sp)==KPROD) | (kt(den_sp)==KSPH)) ? IHAZD : INVLD ); + } + + if ((ang) && (itype != IMULT)) return(INVLD); + + switch(itype) + { case IMULT: + if (ker(den_sp)==WGAUS) return(deg(den_sp)==2); + return( iscompact(ker(den_sp)) ? IMULT : INVLD ); + case IPROD: return( ((den_lfd->d==1) | (kt(den_sp)==KPROD)) ? IPROD : INVLD ); + case IMLIN: return( ((kt(den_sp)==KSPH) && (!lset) && + (deg(den_sp)<=1)) ? IMLIN : INVLD ); + } + + return(INVLD); +} + +extern double lf_tol; + +int densinit(lfd,des,sp) +lfdata *lfd; +design *des; +smpar *sp; +{ int p, i, ii, j, nnz, rnz, ang, lset, status; + double w, *cf; + + den_lfd = lfd; + den_des = des; + den_sp = sp; + cf = des->cf; + + lf_tol = (link(sp)==LLOG) ? 1.0e-6 : 0.0; + + p = des->p; + ff = des->xtwx.wk; + cf[0] = NOSLN; + for (i=1; i<p; i++) cf[i] = 0.0; + + if (!inre(des->xev,lfd->xl,lfd->d)) return(LF_XOOR); + + status = setintlimits(lfd,des->xev,des->h,&ang,&lset); + if (status != LF_OK) return(status); + + switch(selectintmeth(de_itype,lset,ang)) + { case IMULT: des->itype = multint; break; + case IPROD: des->itype = prodint; break; + case IMLIN: des->itype = mlinint; break; + case IHAZD: des->itype = hazint; break; + case INVLD: LERR(("Invalid integration method %d",de_itype)); + break; + case IDEFA: LERR(("No integration type available for this model")); + break; + default: LERR(("densinit: unknown integral type")); + } + + switch(deg(den_sp)) + { case 0: rnz = 1; break; + case 1: rnz = 1; break; + case 2: rnz = lfd->d+1; break; + case 3: rnz = lfd->d+2; break; + default: LERR(("densinit: invalid degree %d",deg(den_sp))); + } + if (lf_error) return(LF_ERR); + + setzero(des->ss,p); + nnz = 0; + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + if (!cens(lfd,ii)) + { w = wght(des,ii)*prwt(lfd,ii); + for (j=0; j<p; j++) des->ss[j] += d_xij(des,ii,j)*w; + if (wght(des,ii)>0.00001) nnz++; + } } + + if (fam(den_sp)==THAZ) haz_init(lfd,des,sp,ilim); +/* this should really only be done once. Not sure how to enforce that, + * esp. when locfit() has been called directly. + */ + if (fam(den_sp)==TDEN) + des->smwt = (lfd->w==NULL) ? lfd->n : vecsum(lfd->w,lfd->n); + + if (lf_debug>2) + { mut_printf(" LHS: "); + for (i=0; i<p; i++) mut_printf(" %8.5f",des->ss[i]); + mut_printf("\n"); + } + + switch(link(den_sp)) + { case LIDENT: + cf[0] = 0.0; + return(LF_OK); + case LLOG: + if (nnz<rnz) { cf[0] = -1000; return(LF_DNOP); } + cf[0] = 0.0; + return(LF_OK); + default: + LERR(("unknown link in densinit")); + return(LF_ERR); + } +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +int bino_vallink(link) +int link; +{ return((link==LLOGIT) | (link==LIDENT) | (link==LASIN)); +} + +int bino_fam(y,p,th,link,res,cens,w) +double y, p, th, *res, w; +int link, cens; +{ double wp; + if (link==LINIT) + { if (y<0) y = 0; + if (y>w) y = w; + res[ZDLL] = y; + return(LF_OK); + } + wp = w*p; + if (link==LIDENT) + { if ((p<=0) && (y>0)) return(LF_BADP); + if ((p>=1) && (y<w)) return(LF_BADP); + res[ZLIK] = res[ZDLL] = res[ZDDLL] = 0.0; + if (y>0) + { res[ZLIK] += y*log(wp/y); + res[ZDLL] += y/p; + res[ZDDLL]+= y/(p*p); + } + if (y<w) + { res[ZLIK] += (w-y)*log((w-wp)/(w-y)); + res[ZDLL] -= (w-y)/(1-p); + res[ZDDLL]+= (w-y)/SQR(1-p); + } + return(LF_OK); + } + if (link==LLOGIT) + { if ((y<0) | (y>w)) /* goon observation; delete it */ + { res[ZLIK] = res[ZDLL] = res[ZDDLL] = 0.0; + return(LF_OK); + } + res[ZLIK] = (th<0) ? th*y-w*log(1+exp(th)) : th*(y-w)-w*log(1+exp(-th)); + if (y>0) res[ZLIK] -= y*log(y/w); + if (y<w) res[ZLIK] -= (w-y)*log(1-y/w); + res[ZDLL] = (y-wp); + res[ZDDLL]= wp*(1-p); + return(LF_OK); + } + if (link==LASIN) + { if ((p<=0) && (y>0)) return(LF_BADP); + if ((p>=1) && (y<w)) return(LF_BADP); + if ((th<0) | (th>PI/2)) return(LF_BADP); + res[ZDLL] = res[ZDDLL] = res[ZLIK] = 0; + if (y>0) + { res[ZDLL] += 2*y*sqrt((1-p)/p); + res[ZLIK] += y*log(wp/y); + } + if (y<w) + { res[ZDLL] -= 2*(w-y)*sqrt(p/(1-p)); + res[ZLIK] += (w-y)*log((w-wp)/(w-y)); + } + res[ZDDLL] = 4*w; + return(LF_OK); + } + LERR(("link %d invalid for binomial family",link)); + return(LF_LNK); +} + +int bino_check(sp,des,lfd) +smpar *sp; +design *des; +lfdata *lfd; +{ int i, ii; + double t0, t1; + + if (fabs(des->cf[0])>700) return(LF_OOB); + + /* check for separation. + * this won't detect separation if there's boundary points with + * both 0 and 1 responses. + */ + t0 = -1e100; t1 = 1e100; + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + if ((resp(lfd,ii)<prwt(lfd,ii)) && (fitv(des,ii) > t0)) t0 = fitv(des,ii); + if ((resp(lfd,ii)>0) && (fitv(des,ii) < t1)) t1 = fitv(des,ii); + if (t1 <= t0) return(LF_OK); + } + mut_printf("separated %8.5f %8.5f\n",t0,t1); + return(LF_NSLN); +} + +void setfbino(fam) +family *fam; +{ fam->deflink = LLOGIT; + fam->canlink = LLOGIT; + fam->vallink = bino_vallink; + fam->family = bino_fam; + fam->pcheck = bino_check; +} + +int rbin_vallink(link) +int link; +{ return(link==LLOGIT); +} + +int rbin_fam(y,p,th,link,res,cens,w) +double y, p, th, *res, w; +int link, cens; +{ double s2y; + if (link==LINIT) + { res[ZDLL] = y; + return(LF_OK); + } + if ((y<0) | (y>w)) /* goon observation; delete it */ + { res[ZLIK] = res[ZDLL] = res[ZDDLL] = 0.0; + return(LF_OK); + } + res[ZLIK] = (th<0) ? th*y-w*log(1+exp(th)) : th*(y-w)-w*log(1+exp(-th)); + if (y>0) res[ZLIK] -= y*log(y/w); + if (y<w) res[ZLIK] -= (w-y)*log(1-y/w); + res[ZDLL] = (y-w*p); + res[ZDDLL]= w*p*(1-p); + if (-res[ZLIK]>HUBERC*HUBERC/2.0) + { s2y = sqrt(-2*res[ZLIK]); + res[ZLIK] = HUBERC*(HUBERC/2.0-s2y); + res[ZDLL] *= HUBERC/s2y; + res[ZDDLL] = HUBERC/s2y*(res[ZDDLL]-1/(s2y*s2y)*w*p*(1-p)); + } + return(LF_OK); +} + +void setfrbino(fam) +family *fam; +{ fam->deflink = LLOGIT; + fam->canlink = LLOGIT; + fam->vallink = rbin_vallink; + fam->family = rbin_fam; + fam->pcheck = bino_check; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +int circ_vallink(link) +int link; +{ return(link==LIDENT); +} + +int circ_fam(y,mean,th,link,res,cens,w) +double y, mean, th, *res, w; +int link, cens; +{ if (link==LINIT) + { res[ZDLL] = w*sin(y); + res[ZLIK] = w*cos(y); + return(LF_OK); + } + res[ZDLL] = w*sin(y-mean); + res[ZDDLL]= w*cos(y-mean); + res[ZLIK] = res[ZDDLL]-w; + return(LF_OK); +} + +extern double lf_tol; +int circ_init(lfd,des,sp) +lfdata *lfd; +design *des; +smpar *sp; +{ int i, ii; + double s0, s1; + s0 = s1 = 0.0; + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + s0 += wght(des,ii)*prwt(lfd,ii)*sin(resp(lfd,ii)-base(lfd,ii)); + s1 += wght(des,ii)*prwt(lfd,ii)*cos(resp(lfd,ii)-base(lfd,ii)); + } + des->cf[0] = atan2(s0,s1); + for (i=1; i<des->p; i++) des->cf[i] = 0.0; + lf_tol = 1.0e-6; + return(LF_OK); +} + + +void setfcirc(fam) +family *fam; +{ fam->deflink = LIDENT; + fam->canlink = LIDENT; + fam->vallink = circ_vallink; + fam->family = circ_fam; + fam->initial = circ_init; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +int dens_vallink(link) +int link; +{ return((link==LIDENT) | (link==LLOG)); +} + +int dens_fam(y,mean,th,link,res,cens,w) +double y, mean, th, *res, w; +int link, cens; +{ if (cens) + res[ZLIK] = res[ZDLL] = res[ZDDLL] = 0.0; + else + { res[ZLIK] = w*th; + res[ZDLL] = res[ZDDLL] = w; + } + return(LF_OK); +} + +void setfdensity(fam) +family *fam; +{ fam->deflink = LLOG; + fam->canlink = LLOG; + fam->vallink = dens_vallink; + fam->family = dens_fam; + fam->initial = densinit; + fam->like = likeden; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +int gamma_vallink(link) +int link; +{ return((link==LIDENT) | (link==LLOG) | (link==LINVER)); +} + +int gamma_fam(y,mean,th,link,res,cens,w) +double y, mean, th, *res, w; +int link, cens; +{ double lb, pt, dg; + if (link==LINIT) + { res[ZDLL] = MAX(y,0.0); + return(LF_OK); + } + res[ZLIK] = res[ZDLL] = res[ZDDLL] = 0.0; + if (w==0.0) return(LF_OK); + if ((mean<=0) & (y>0)) return(LF_BADP); + if (link==LIDENT) lb = 1/th; + if (link==LINVER) lb = th; + if (link==LLOG) lb = mut_exp(-th); + if (cens) + { if (y<=0) return(LF_OK); + pt = 1-igamma(lb*y,w); + dg = dgamma(lb*y,w,1.0,0); + res[ZLIK] = log(pt); + res[ZDLL] = -y*dg/pt; +/* + * res[ZDLL] = -y*dg/pt * dlb/dth. + * res[ZDDLL] = y*dg/pt * (d2lb/dth2 + ((w-1)/lb-y)*(dlb/dth)^2) + * + res[ZDLL]^2. + */ + if (link==LLOG) /* lambda = exp(-theta) */ + { res[ZDLL] *= -lb; + res[ZDDLL] = dg*y*lb*(w-lb*y)/pt + SQR(res[ZDLL]); + return(LF_OK); + } + if (link==LINVER) /* lambda = theta */ + { res[ZDLL] *= 1.0; + res[ZDDLL] = dg*y*((w-1)*mean-y)/pt + SQR(res[ZDLL]); + return(LF_OK); + } + if (link==LIDENT) /* lambda = 1/theta */ + { res[ZDLL] *= -lb*lb; + res[ZDDLL] = dg*y*lb*lb*lb*(1+w-lb*y)/pt + SQR(res[ZDLL]); + return(LF_OK); + } + } + else + { if (y<0) WARN(("Negative Gamma observation")); + if (link==LLOG) + { res[ZLIK] = -lb*y+w*(1-th); + if (y>0) res[ZLIK] += w*log(y/w); + res[ZDLL] = lb*y-w; + res[ZDDLL]= lb*y; + return(LF_OK); + } + if (link==LINVER) + { res[ZLIK] = -lb*y+w-w*log(mean); + if (y>0) res[ZLIK] += w*log(y/w); + res[ZDLL] = -y+w*mean; + res[ZDDLL]= w*mean*mean; + return(LF_OK); + } + if (link==LIDENT) + { res[ZLIK] = -lb*y+w-w*log(mean); + if (y>0) res[ZLIK] += w*log(y/w); + res[ZDLL] = lb*lb*(y-w*mean); + res[ZDDLL]= lb*lb*lb*(2*y-w*mean); + return(LF_OK); + } + } + LERR(("link %d invalid for Gamma family",link)); + return(LF_LNK); +} + +void setfgamma(fam) +family *fam; +{ fam->deflink = LLOG; + fam->canlink = LINVER; + fam->vallink = gamma_vallink; + fam->family = gamma_fam; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +int gaus_vallink(link) +int link; +{ return((link==LIDENT) | (link==LLOG) | (link==LLOGIT)); +} + +int gaus_fam(y,mean,th,link,res,cens,w) +double y, mean, th, *res, w; +int link, cens; +{ double z, pz, dp; + if (link==LINIT) + { res[ZDLL] = w*y; + return(LF_OK); + } + z = y-mean; + if (cens) + { if (link!=LIDENT) + { LERR(("Link invalid for censored Gaussian family")); + return(LF_LNK); + } + pz = mut_pnorm(-z); + dp = ((z>6) ? ptail(-z) : exp(-z*z/2)/pz)/2.5066283; + res[ZLIK] = w*log(pz); + res[ZDLL] = w*dp; + res[ZDDLL]= w*dp*(dp-z); + return(LF_OK); + } + res[ZLIK] = -w*z*z/2; + switch(link) + { case LIDENT: + res[ZDLL] = w*z; + res[ZDDLL]= w; + break; + case LLOG: + res[ZDLL] = w*z*mean; + res[ZDDLL]= w*mean*mean; + break; + case LLOGIT: + res[ZDLL] = w*z*mean*(1-mean); + res[ZDDLL]= w*mean*mean*(1-mean)*(1-mean); + break; + default: + LERR(("Invalid link for Gaussian family")); + return(LF_LNK); + } + return(LF_OK); +} + +int gaus_check(sp,des,lfd) +smpar *sp; +design *des; +lfdata *lfd; +{ int i, ii; + if (fami(sp)->robust) return(LF_OK); + if (link(sp)==LIDENT) + { for (i=0; i<des->n; i++) + { ii = des->ind[i]; + if (cens(lfd,ii)) return(LF_OK); + } + return(LF_DONE); + } + return(LF_OK); +} + +void setfgauss(fam) +family *fam; +{ fam->deflink = LIDENT; + fam->canlink = LIDENT; + fam->vallink = gaus_vallink; + fam->family = gaus_fam; + fam->pcheck = gaus_check; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +int geom_vallink(link) +int link; +{ return((link==LIDENT) | (link==LLOG)); +} + +int geom_fam(y,mean,th,link,res,cens,w) +double y, mean, th, *res, w; +int link, cens; +{ double p, pt, dp, p1; + if (link==LINIT) + { res[ZDLL] = MAX(y,0.0); + return(LF_OK); + } + p = 1/(1+mean); + if (cens) /* censored observation */ + { if (y<=0) + { res[ZLIK] = res[ZDLL] = res[ZDDLL] = 0; + return(LF_OK); + } + p1 = (link==LIDENT) ? -p*p : -p*(1-p); + pt = 1-ibeta(p,w,y); + dp = dbeta(p,w,y,0)/pt; + res[ZLIK] = log(pt); + res[ZDLL] = -dp*p1; + res[ZDDLL] = dp*dp*p1*p1; + if (link==LIDENT) + res[ZDDLL] += dp*p*p*p*(1+w*(1-p)-p*y)/(1-p); + else + res[ZDDLL] += dp*p*(1-p)*(w*(1-p)-p*y); + return(LF_OK); + } + else + { res[ZLIK] = (y+w)*log((y/w+1)/(mean+1)); + if (y>0) res[ZLIK] += y*log(w*mean/y); + if (link==LLOG) + { res[ZDLL] = (y-w*mean)*p; + res[ZDDLL]= (y+w)*p*(1-p); + return(LF_OK); + } + if (link==LIDENT) + { res[ZDLL] = (y-w*mean)/(mean*(1+mean)); + res[ZDDLL]= w/(mean*(1+mean)); + return(LF_OK); + } + } + LERR(("link %d invalid for geometric family",link)); + return(LF_LNK); +} + +void setfgeom(fam) +family *fam; +{ fam->deflink = LLOG; + fam->canlink = LIDENT; /* this isn't correct. I haven't prog. canon */ + fam->vallink = geom_vallink; + fam->family = geom_fam; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +#define HUBERC 2.0 + +double links_rs; +int inllmix=0; + +/* + * lffamily("name") converts family names into a numeric value. + * typical usage is fam(&lf->sp) = lffamily("gaussian"); + * Note that family can be preceded by q and/or r for quasi, robust. + * + * link(&lf->sp) = lflink("log") does the same for the link function. + */ +#define NFAMILY 18 +static char *famil[NFAMILY] = + { "density", "ate", "hazard", "gaussian", "binomial", + "poisson", "gamma", "geometric", "circular", "obust", "huber", + "weibull", "cauchy","probab", "logistic", "nbinomial", + "vonmises", "quant" }; +static int fvals[NFAMILY] = + { TDEN, TRAT, THAZ, TGAUS, TLOGT, + TPOIS, TGAMM, TGEOM, TCIRC, TROBT, TROBT, + TWEIB, TCAUC, TPROB, TLOGT, TGEOM, TCIRC, TQUANT }; +int lffamily(z) +char *z; +{ int quasi, robu, f; + quasi = robu = 0; + while ((z[0]=='q') | (z[0]=='r')) + { quasi |= (z[0]=='q'); + robu |= (z[0]=='r'); + z++; + } + z[0] = tolower(z[0]); + f = pmatch(z,famil,fvals,NFAMILY,-1); + if ((z[0]=='o') | (z[0]=='a')) robu = 0; + if (f==-1) + { WARN(("unknown family %s",z)); + f = TGAUS; + } + if (quasi) f += 64; + if (robu) f += 128; + return(f); +} + +#define NLINKS 8 +static char *ltype[NLINKS] = { "default", "canonical", "identity", "log", + "logi", "inverse", "sqrt", "arcsin" }; +static int lvals[NLINKS] = { LDEFAU, LCANON, LIDENT, LLOG, + LLOGIT, LINVER, LSQRT, LASIN }; +int lflink(char *z) +{ int f; + if (z==NULL) return(LDEFAU); + z[0] = tolower(z[0]); + f = pmatch(z, ltype, lvals, NLINKS, -1); + if (f==-1) + { WARN(("unknown link %s",z)); + f = LDEFAU; + } + return(f); +} + +int defaultlink(link,fam) +int link; +family *fam; +{ if (link==LDEFAU) return(fam->deflink); + if (link==LCANON) return(fam->canlink); + return(link); +} + +/* +void robustify(res,rs) +double *res, rs; +{ double sc, z; + sc = rs*HUBERC; + if (res[ZLIK] > -sc*sc/2) return; + z = sqrt(-2*res[ZLIK]); + res[ZDDLL]= -sc*res[ZDLL]*res[ZDLL]/(z*z*z)+sc*res[ZDDLL]/z; + res[ZDLL]*= sc/z; + res[ZLIK] = sc*sc/2-sc*z; +} +*/ +void robustify(res,rs) +double *res, rs; +{ double sc, z; + sc = rs*HUBERC; + if (res[ZLIK] > -sc*sc/2) + { res[ZLIK] /= sc*sc; + res[ZDLL] /= sc*sc; + res[ZDDLL] /= sc*sc; + return; + } + z = sqrt(-2*res[ZLIK]); + res[ZDDLL]= (-sc*res[ZDLL]*res[ZDLL]/(z*z*z)+sc*res[ZDDLL]/z)/(sc*sc); + res[ZDLL]*= 1.0/(z*sc); + res[ZLIK] = 0.5-z/sc; +} + +double lf_link(y,lin) +double y; +int lin; +{ switch(lin) + { case LIDENT: return(y); + case LLOG: return(log(y)); + case LLOGIT: return(logit(y)); + case LINVER: return(1/y); + case LSQRT: return(sqrt(fabs(y))); + case LASIN: return(asin(sqrt(y))); + } + LERR(("link: unknown link %d",lin)); + return(0.0); +} + +double invlink(th,lin) +double th; +int lin; +{ switch(lin) + { case LIDENT: return(th); + case LLOG: return(mut_exp(th)); + case LLOGIT: return(expit(th)); + case LINVER: return(1/th); + case LSQRT: return(th*fabs(th)); + case LASIN: return(sin(th)*sin(th)); + case LINIT: return(0.0); + } + LERR(("invlink: unknown link %d",lin)); + return(0.0); +} + +/* the link and various related functions */ +int links(th,y,fam,link,res,c,w,rs) +double th, y, *res, w, rs; +int link, c; +family *fam; +{ double mean; + int st; + + mean = res[ZMEAN] = invlink(th,link); + if (lf_error) return(LF_LNK); + links_rs = rs; +/* mut_printf("links: rs %8.5f\n",rs); */ + + st = fam->family(y,mean,th,link,res,c,w); + + if (st!=LF_OK) return(st); + if (link==LINIT) return(st); + if (isrobust(fam)) robustify(res,rs); + return(st); +} + +/* + stdlinks is a version of links when family, link, response e.t.c + all come from the standard places. +*/ +int stdlinks(res,lfd,sp,i,th,rs) +lfdata *lfd; +smpar *sp; +double th, rs, *res; +int i; +{ + return(links(th,resp(lfd,i),fami(sp),link(sp),res,cens(lfd,i),prwt(lfd,i),rs)); +} + +/* + * functions used in variance, skewness, kurtosis calculations + * in scb corrections. + */ + +double b2(th,tg,w) +double th, w; +int tg; +{ double y; + switch(tg&63) + { case TGAUS: return(w); + case TPOIS: return(w*mut_exp(th)); + case TLOGT: + y = expit(th); + return(w*y*(1-y)); + } + LERR(("b2: invalid family %d",tg)); + return(0.0); +} + +double b3(th,tg,w) +double th, w; +int tg; +{ double y; + switch(tg&63) + { case TGAUS: return(0.0); + case TPOIS: return(w*mut_exp(th)); + case TLOGT: + y = expit(th); + return(w*y*(1-y)*(1-2*y)); + } + LERR(("b3: invalid family %d",tg)); + return(0.0); +} + +double b4(th,tg,w) +double th, w; +int tg; +{ double y; + switch(tg&63) + { case TGAUS: return(0.0); + case TPOIS: return(w*mut_exp(th)); + case TLOGT: + y = expit(th); y = y*(1-y); + return(w*y*(1-6*y)); + } + LERR(("b4: invalid family %d",tg)); + return(0.0); +} + +int def_check(sp,des,lfd) +smpar *sp; +design *des; +lfdata *lfd; +{ switch(link(sp)) + { case LLOG: if (des->cf[0]>700) return(LF_OOB); + break; + } + return(LF_OK); +} +extern void setfdensity(), setfgauss(), setfbino(), setfpoisson(); +extern void setfgamma(), setfgeom(), setfcirc(), setfweibull(); +extern void setfrbino(), setfrobust(), setfcauchy(), setfquant(); + +void setfamily(sp) +smpar *sp; +{ int tg, lnk; + family *f; + + tg = fam(sp); + f = fami(sp); + f->quasi = tg&64; + f->robust = tg&128; + f->initial = reginit; + f->like = likereg; + f->pcheck = def_check; + + switch(tg&63) + { case TDEN: + case THAZ: + case TRAT: setfdensity(f); break; + case TGAUS: setfgauss(f); break; + case TLOGT: setfbino(f); break; + case TRBIN: setfrbino(f); break; + case TPROB: + case TPOIS: setfpoisson(f); break; + case TGAMM: setfgamma(f); break; + case TGEOM: setfgeom(f); break; + case TWEIB: setfweibull(f); + case TCIRC: setfcirc(f); break; + case TROBT: setfrobust(f); break; + case TCAUC: setfcauchy(f); break; + case TQUANT: setfquant(f); break; + default: LERR(("setfamily: unknown family %d",tg&63)); + return; + } + + lnk = defaultlink(link(sp),f); + if (!f->vallink(lnk)) + { WARN(("setfamily: invalid link %d - revert to default",link(sp))); + link(sp) = f->deflink; + } + else + link(sp) = lnk; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +int pois_vallink(link) +int link; +{ return((link==LLOG) | (link==LIDENT) | (link==LSQRT)); +} + +int pois_fam(y,mean,th,link,res,cens,w) +double y, mean, th, *res, w; +int link, cens; +{ double wmu, pt, dp; + if (link==LINIT) + { res[ZDLL] = MAX(y,0.0); + return(LF_OK); + } + wmu = w*mean; + if (inllmix) y = w*y; + if (cens) + { if (y<=0) + { res[ZLIK] = res[ZDLL] = res[ZDDLL] = 0.0; + return(LF_OK); + } + pt = igamma(wmu,y); + dp = dgamma(wmu,y,1.0,0)/pt; + res[ZLIK] = log(pt); +/* + * res[ZDLL] = dp * w*dmu/dth + * res[ZDDLL]= -dp*(w*d2mu/dth2 + (y-1)/mu*(dmu/dth)^2) + res[ZDLL]^2 + */ + if (link==LLOG) + { res[ZDLL] = dp*wmu; + res[ZDDLL]= -dp*wmu*(y-wmu) + SQR(res[ZDLL]); + return(LF_OK); + } + if (link==LIDENT) + { res[ZDLL] = dp*w; + res[ZDDLL]= -dp*(y-1-wmu)*w/mean + SQR(res[ZDLL]); + return(LF_OK); + } + if (link==LSQRT) + { res[ZDLL] = dp*2*w*th; + res[ZDDLL]= -dp*w*(4*y-2-4*wmu) + SQR(res[ZDLL]); + return(LF_OK); + } } + if (link==LLOG) + { if (y<0) /* goon observation - delete it */ + { res[ZLIK] = res[ZDLL] = res[ZDDLL] = 0; + return(LF_OK); + } + res[ZLIK] = res[ZDLL] = y-wmu; + if (y>0) res[ZLIK] += y*(th-log(y/w)); + res[ZDDLL] = wmu; + return(LF_OK); + } + if (link==LIDENT) + { if ((mean<=0) && (y>0)) return(LF_BADP); + res[ZLIK] = y-wmu; + res[ZDLL] = -w; + res[ZDDLL] = 0; + if (y>0) + { res[ZLIK] += y*log(wmu/y); + res[ZDLL] += y/mean; + res[ZDDLL]= y/(mean*mean); + } + return(LF_OK); + } + if (link==LSQRT) + { if ((mean<=0) && (y>0)) return(LF_BADP); + res[ZLIK] = y-wmu; + res[ZDLL] = -2*w*th; + res[ZDDLL]= 2*w; + if (y>0) + { res[ZLIK] += y*log(wmu/y); + res[ZDLL] += 2*y/th; + res[ZDDLL]+= 2*y/mean; + } + return(LF_OK); + } + LERR(("link %d invalid for Poisson family",link)); + return(LF_LNK); +} + +void setfpoisson(fam) +family *fam; +{ fam->deflink = LLOG; + fam->canlink = LLOG; + fam->vallink = pois_vallink; + fam->family = pois_fam; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +#define QTOL 1.0e-10 +extern int lf_status; +static double q0; + +int quant_vallink(int link) { return(1); } + +int quant_fam(y,mean,th,link,res,cens,w) +double y, mean, th, *res, w; +int link, cens; +{ double z, p; + if (link==LINIT) + { res[ZDLL] = w*y; + return(LF_OK); + } +p = 0.5; /* should be pen(sp) */ + z = y-mean; + res[ZLIK] = (z<0) ? (w*z/p) : (-w*z/(1-p)); + res[ZDLL] = (z<0) ? -w/p : w/(1-p); + res[ZDDLL]= w/(p*(1-p)); + return(LF_OK); +} + +int quant_check(sp,des,lfd) +smpar *sp; +design *des; +lfdata *lfd; +{ return(LF_DONE); +} + +void setfquant(fam) +family *fam; +{ fam->deflink = LIDENT; + fam->canlink = LIDENT; + fam->vallink = quant_vallink; + fam->family = quant_fam; + fam->pcheck = quant_check; +} + +/* + * cycling rule for choosing among ties. + */ +int tiecycle(ind,i0,i1,oi) +int *ind, i0, i1, oi; +{ int i, ii, im; + im = ind[i0]; + for (i=i0+1; i<=i1; i++) + { ii = ind[i]; + if (im<=oi) + { if ((ii<im) | (ii>oi)) im = ii; + } + else + { if ((ii<im) & (ii>oi)) im = ii; + } + } + return(im); +} + +/* + * move coefficient vector cf, as far as possible, in direction dc. + */ +int movecoef(lfd,des,p,cf,dc,oi) +lfdata *lfd; +design *des; +double p, *cf, *dc; +int oi; +{ int i, ii, im, i0, i1, j; + double *lb, *el, e, sp, sn, sw, sum1, sum2, tol1; + + lb = des->th; + el = des->res; + sum1 = sum2 = 0.0; + + sp = sn = sw = 0.0; + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + lb[ii] = innerprod(dc,d_xi(des,ii),des->p); + e = resp(lfd,ii) - innerprod(cf,d_xi(des,ii),des->p); + el[ii] = (fabs(lb[ii])<QTOL) ? 1e100 : e/lb[ii]; + if (lb[ii]>0) + sp += prwt(lfd,ii)*wght(des,ii)*lb[ii]; + else + sn -= prwt(lfd,ii)*wght(des,ii)*lb[ii]; + sw += prwt(lfd,ii)*wght(des,ii); + } +printf("sp %8.5f sn %8.5f\n",sn,sp); +/* if sn, sp are both zero, should return an LF_PF. + * but within numerical tolerance? what does it mean? + */ + if (sn+sp <= QTOL*q0) { lf_status = LF_PF; return(0); } + + sum1 = sp/(1-p) + sn/p; + tol1 = QTOL*(sp+sn); + mut_order(el,des->ind,0,des->n-1); + + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + sum2 += prwt(lfd,ii)*wght(des,ii)*((lb[ii]>0) ? lb[ii]/p : -lb[ii]/(1-p) ); + sum1 -= prwt(lfd,ii)*wght(des,ii)*((lb[ii]>0) ? lb[ii]/(1-p) : -lb[ii]/p ); + if (sum1<=sum2+tol1) + { +/* determine the range of ties [i0,i1] + * el[ind[i0..i1]] = el[ind[i]]. + * if sum1==sum2, el[ind[i+1]]..el[ind[i1]]] = el[ind[i1]], else i1 = i. + */ + i0 = i1 = i; + while ((i0>0) && (el[des->ind[i0-1]]==el[ii])) i0--; + while ((i1<des->n-1) && (el[des->ind[i1+1]]==el[ii])) i1++; + if (sum1>=sum2-tol1) + while ((i1<des->n-1) && (el[des->ind[i1+1]]==el[des->ind[i+1]])) i1++; + + if (i0<i1) ii = tiecycle(des->ind,i0,i1,oi); + for (j=0; j<des->p; j++) cf[j] += el[ii]*dc[j]; + return(ii); + } + } +mut_printf("Big finddlt problem.\n"); +ii = des->ind[des->n-1]; +for (j=0; j<des->p; j++) cf[j] += el[ii]*dc[j]; +return(ii); +} + +/* + * special version of movecoef for min/max. + */ +int movemin(lfd,des,f,cf,dc,oi) +design *des; +lfdata *lfd; +double *cf, *dc, f; +int oi; +{ int i, ii, im, p, s, ssum; + double *lb, sum, lb0, lb1, z0, z1; + + lb = des->th; + s = (f<=0.0) ? 1 : -1; + +/* first, determine whether move should be in positive or negative direction */ + p = des->p; + sum = 0; + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + lb[ii] = innerprod(dc,d_xi(des,ii),des->p); + sum += prwt(lfd,ii)*wght(des,ii)*lb[ii]; + } + if (fabs(sum) <= QTOL*q0) + { lf_status = LF_PF; + return(0); + } + ssum = (sum<=0.0) ? -1 : 1; + if (ssum != s) + for (i=0; i<p; i++) dc[i] = -dc[i]; + +/* now, move positively. How far can we move? */ + lb0 = 1.0e100; im = oi; + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + lb[ii] = innerprod(dc,d_xi(des,ii),des->p); /* must recompute - signs! */ + if (s*lb[ii]>QTOL) /* should have scale-free tolerance here */ + { z0 = innerprod(cf,d_xi(des,ii),p); + lb1 = (resp(lfd,ii) - z0)/lb[ii]; + if (lb1<lb0) + { if (fabs(lb1-lb0)<QTOL) /* cycle */ + { if (im<=oi) + { if ((ii>oi) | (ii<im)) im = ii; } + else + { if ((ii>oi) & (ii<im)) im = ii; } + } + else + { im = ii; lb0 = lb1; } + } + } + } + + for (i=0; i<p; i++) cf[i] = cf[i]+lb0*dc[i]; + if (im==-1) lf_status = LF_PF; + return(im); +} + +double qll(lfd,spr,des,cf) +lfdata *lfd; +smpar *spr; +design *des; +double *cf; +{ int i, ii; + double th, sp, sn, p, e; + + p = pen(spr); + sp = sn = 0.0; + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + th = innerprod(d_xi(des,ii),cf,des->p); + e = resp(lfd,ii)-th; + if (e<0) sn -= prwt(lfd,ii)*wght(des,ii)*e; + if (e>0) sp += prwt(lfd,ii)*wght(des,ii)*e; + } + if (p<=0.0) return((sn<QTOL) ? -sp : -1e300); + if (p>=1.0) return((sp<QTOL) ? -sn : -1e300); + return(-sp/(1-p)-sn/p); +} + +/* + * running quantile smoother. + */ +void lfquantile(lfd,sp,des,maxit) +lfdata *lfd; +smpar *sp; +design *des; +int maxit; +{ int i, ii, im, j, k, p, *ci, (*mover)(); + double *cf, *db, *dc, *cm, f, q1, q2, l0; + +printf("in lfquantile\n"); + f = pen(sp); + p = des->p; + cf = des->cf; + dc = des->oc; + db = des->ss; + setzero(cf,p); + setzero(dc,p); + cm = des->V; + setzero(cm,p*p); + ci = (int *)des->fix; + + q1 = -qll(lfd,sp,des,cf); + if (q1==0.0) { lf_status = LF_PF; return; } + for (i=0; i<p; i++) cm[i*(p+1)] = 1; + mover = movecoef; + if ((f<=0.0) | (f>=1.0)) mover = movemin; + + dc[0] = 1.0; + im = mover(lfd,des,f,cf,dc,-1); + if (lf_status != LF_OK) return; + ci[0] = im; +printf("init const %2d\n",ci[0]); + q0 = -qll(lfd,sp,des,cf); + if (q0<QTOL*q1) { lf_status = LF_PF; return; } + +printf("loop 0\n"); fflush(stdout); + for (i=1; i<p; i++) + { +printf("i %2d\n",i); + memcpy(&cm[(i-1)*p],d_xi(des,im),p*sizeof(double)); + setzero(db,p); + db[i] = 1.0; + resproj(db,cm,dc,p,i); +printf("call mover\n"); fflush(stdout); + im = mover(lfd,des,f,cf,dc,-1); + if (lf_status != LF_OK) return; +printf("mover %2d\n",im); fflush(stdout); + ci[i] = im; + } +printf("call qll\n"); fflush(stdout); + q1 = qll(lfd,sp,des,cf); + +printf("loop 1 %d %d %d %d\n",ci[0],ci[1],ci[2],ci[3]); fflush(stdout); + for (k=0; k<maxit; k++) + { for (i=0; i<p; i++) + { for (j=0; j<p; j++) + if (j!=i) memcpy(&cm[(j-(j>i))*p],d_xi(des,ci[j]),p*sizeof(double)); + memcpy(db,d_xi(des,ci[i]),p*sizeof(double)); + resproj(db,cm,dc,p,p-1); +printf("call mover\n"); fflush(stdout); + im = mover(lfd,des,f,cf,dc,ci[i]); + if (lf_status != LF_OK) return; +printf("mover %2d\n",im); fflush(stdout); + ci[i] = im; + } + q2 = qll(lfd,sp,des,cf); +/* + * convergence: require no change -- reasonable, since discrete? + * remember we're maximizing, and q's are negative. + */ + if (q2 <= q1) return; + q1 = q2; + } +printf("loop 2\n"); + mut_printf("Warning: lfquantile not converged.\n"); +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +extern double links_rs; + +int robust_vallink(link) +int link; +{ return(link==LIDENT); +} + +int robust_fam(y,mean,th,link,res,cens,w) +double y, mean, th, *res, w; +int link, cens; +{ double z, sw; + if (link==LINIT) + { res[ZDLL] = w*y; + return(LF_OK); + } + sw = (w==1.0) ? 1.0 : sqrt(w); /* don't want unnecess. sqrt! */ + z = sw*(y-mean)/links_rs; + res[ZLIK] = (fabs(z)<HUBERC) ? -z*z/2 : HUBERC*(HUBERC/2.0-fabs(z)); + if (z< -HUBERC) + { res[ZDLL] = -sw*HUBERC/links_rs; + res[ZDDLL]= 0.0; + return(LF_OK); + } + if (z> HUBERC) + { res[ZDLL] = sw*HUBERC/links_rs; + res[ZDDLL]= 0.0; + return(LF_OK); + } + res[ZDLL] = sw*z/links_rs; + res[ZDDLL] = w/(links_rs*links_rs); + return(LF_OK); +} + +int cauchy_fam(y,p,th,link,res,cens,w) +double y, p, th, *res, w; +int link, cens; +{ double z; + if (link!=LIDENT) + { LERR(("Invalid link in famcauc")); + return(LF_LNK); + } + z = w*(y-th)/links_rs; + res[ZLIK] = -log(1+z*z); + res[ZDLL] = 2*w*z/(links_rs*(1+z*z)); + res[ZDDLL] = 2*w*w*(1-z*z)/(links_rs*links_rs*(1+z*z)*(1+z*z)); + return(LF_OK); +} + +extern double lf_tol; +int robust_init(lfd,des,sp) +lfdata *lfd; +design *des; +smpar *sp; +{ int i; + for (i=0; i<des->n; i++) + des->res[i] = resp(lfd,(int)des->ind[i]) - base(lfd,(int)des->ind[i]); + des->cf[0] = median(des->res,des->n); + for (i=1; i<des->p; i++) des->cf[i] = 0.0; + lf_tol = 1.0e-6; + return(LF_OK); +} + +void setfrobust(fam) +family *fam; +{ fam->deflink = LIDENT; + fam->canlink = LIDENT; + fam->vallink = robust_vallink; + fam->family = robust_fam; + fam->initial = robust_init; + fam->robust = 0; +} + +void setfcauchy(fam) +family *fam; +{ fam->deflink = LIDENT; + fam->canlink = LIDENT; + fam->vallink = robust_vallink; + fam->family = cauchy_fam; + fam->initial = robust_init; + fam->robust = 0; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +int weibull_vallink(link) +int link; +{ return((link==LIDENT) | (link==LLOG) | (link==LLOGIT)); +} + +int weibull_fam(y,mean,th,link,res,cens,w) +double y, mean, th, *res, w; +int link, cens; +{ double yy; + yy = pow(y,w); + if (link==LINIT) + { res[ZDLL] = MAX(yy,0.0); + return(LF_OK); + } + if (cens) + { res[ZLIK] = -yy/mean; + res[ZDLL] = res[ZDDLL] = yy/mean; + return(LF_OK); + } + res[ZLIK] = 1-yy/mean-th; + if (yy>0) res[ZLIK] += log(w*yy); + res[ZDLL] = -1+yy/mean; + res[ZDDLL]= yy/mean; + return(LF_OK); +} + +void setfweibull(fam) +family *fam; +{ fam->deflink = LLOG; + fam->canlink = LLOG; + fam->vallink = weibull_vallink; + fam->family = weibull_fam; + fam->robust = 0; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + Functions implementing the adaptive bandwidth selection. + Will make the final call to nbhd() to set smoothing weights + for selected bandwidth, But will **not** make the + final call to locfit(). +*/ + +#include "locf.h" + +static double hmin; + +#define NACRI 5 +static char *atype[NACRI] = { "none", "cp", "ici", "mindex", "ok" }; +static int avals[NACRI] = { ANONE, ACP, AKAT, AMDI, AOK }; +int lfacri(char *z) +{ return(pmatch(z, atype, avals, NACRI, ANONE)); +} + +double adcri(lk,t0,t2,pen) +double lk, t0, t2, pen; +{ double y; +/* return(-2*lk/(t0*exp(pen*log(1-t2/t0)))); */ + /* return((-2*lk+pen*t2)/t0); */ + y = (MAX(-2*lk,t0-t2)+pen*t2)/t0; + return(y); +} + +double mmse(lfd,sp,dv,des) +lfdata *lfd; +smpar *sp; +deriv *dv; +design *des; +{ int i, ii, j, p, p1; + double sv, sb, *l, dp; + + l = des->wd; + wdiag(lfd, sp, des,l,dv,0,1,0); + sv = sb = 0; + p = npar(sp); + for (i=0; i<des->n; i++) + { sv += l[i]*l[i]; + ii = des->ind[i]; + dp = dist(des,ii); + for (j=0; j<deg(sp); j++) dp *= dist(des,ii); + sb += fabs(l[i])*dp; + } + p1 = factorial(deg(sp)+1); +printf("%8.5f sv %8.5f sb %8.5f %8.5f\n",des->h,sv,sb,sv+sb*sb*pen(sp)*pen(sp)/(p1*p1)); + return(sv+sb*sb*pen(sp)*pen(sp)/(p1*p1)); +} + +static double mcp, clo, cup; + +/* + Initial bandwidth will be (by default) + k-nearest neighbors for k small, just large enough to + get defined estimate (unless user provided nonzero nn or fix-h components) +*/ + +int ainitband(lfd,sp,dv,des) +lfdata *lfd; +smpar *sp; +deriv *dv; +design *des; +{ int lf_status, p, z, cri, noit, redo; + double ho, t[6]; + + if (lf_debug >= 2) mut_printf("ainitband:\n"); + p = des->p; + cri = acri(sp); + noit = (cri!=AOK); + z = (int)(lfd->n*nn(sp)); + if ((noit) && (z<p+2)) z = p+2; + redo = 0; ho = -1; + do + { + nbhd(lfd,des,z,redo,sp); + if (z<des->n) z = des->n; + if (des->h>ho) lf_status = locfit(lfd,des,sp,noit,0,0); + z++; + redo = 1; + } while ((z<=lfd->n) && ((des->h==0)||(lf_status!=LF_OK))); + hmin = des->h; + + switch(cri) + { case ACP: + local_df(lfd,sp,des,t); + mcp = adcri(des->llk,t[0],t[2],pen(sp)); + return(lf_status); + case AKAT: + local_df(lfd,sp,des,t); + clo = des->cf[0]-pen(sp)*t[5]; + cup = des->cf[0]+pen(sp)*t[5]; + return(lf_status); + case AMDI: + mcp = mmse(lfd,sp,dv,des); + return(lf_status); + case AOK: return(lf_status); + } + LERR(("aband1: unknown criterion")); + return(LF_ERR); +} + +/* + aband2 increases the initial bandwidth until lack of fit results, + or the fit is close to a global fit. Increase h by 1+0.3/d at + each iteration. +*/ + +double aband2(lfd,sp,dv,des,h0) +lfdata *lfd; +smpar *sp; +deriv *dv; +design *des; +double h0; +{ double t[6], h1, nu1, cp, ncp, tlo, tup; + int d, inc, n, p, done; + + if (lf_debug >= 2) mut_printf("aband2:\n"); + d = lfd->d; n = lfd->n; p = npar(sp); + h1 = des->h = h0; + done = 0; nu1 = 0.0; + inc = 0; ncp = 0.0; + while ((!done) & (nu1<(n-p)*0.95)) + { fixh(sp) = (1+0.3/d)*des->h; + nbhd(lfd,des,0,1,sp); + if (locfit(lfd,des,sp,1,0,0) > 0) WARN(("aband2: failed fit")); + local_df(lfd,sp,des,t); + nu1 = t[0]-t[2]; /* tr(A) */ + switch(acri(sp)) + { case AKAT: + tlo = des->cf[0]-pen(sp)*t[5]; + tup = des->cf[0]+pen(sp)*t[5]; +/* mut_printf("h %8.5f tlo %8.5f tup %8.5f\n",des->h,tlo,tup); */ + done = ((tlo>cup) | (tup<clo)); + if (!done) + { clo = MAX(clo,tlo); + cup = MIN(cup,tup); + h1 = des->h; + } + break; + case ACP: + cp = adcri(des->llk,t[0],t[2],pen(sp)); +/* mut_printf("h %8.5f lk %8.5f t0 %8.5f t2 %8.5f cp %8.5f\n",des->h,des->llk,t[0],t[2],cp); */ + if (cp<mcp) { mcp = cp; h1 = des->h; } + if (cp>=ncp) inc++; else inc = 0; + ncp = cp; + done = (inc>=10) | ((inc>=3) & ((t[0]-t[2])>=10) & (cp>1.5*mcp)); + break; + case AMDI: + cp = mmse(lfd,sp,dv,des); + if (cp<mcp) { mcp = cp; h1 = des->h; } + if (cp>ncp) inc++; else inc = 0; + ncp = cp; + done = (inc>=3); + break; + } + } + return(h1); +} + +/* + aband3 does a finer search around best h so far. Try + h*(1-0.2/d), h/(1-0.1/d), h*(1+0.1/d), h*(1+0.2/d) +*/ +double aband3(lfd,sp,dv,des,h0) +lfdata *lfd; +smpar *sp; +deriv *dv; +design *des; +double h0; +{ double t[6], h1, cp, tlo, tup; + int i, i0, d, n; + + if (lf_debug >= 2) mut_printf("aband3:\n"); + d = lfd->d; n = lfd->n; + h1 = h0; + i0 = (acri(sp)==AKAT) ? 1 : -2; + if (h0==hmin) i0 = 1; + + for (i=i0; i<=2; i++) + { if (i==0) i++; + fixh(sp) = h0*(1+0.1*i/d); + nbhd(lfd,des,0,1,sp); + if (locfit(lfd,des,sp,1,0,0) > 0) WARN(("aband3: failed fit")); + local_df(lfd,sp,des,t); + switch (acri(sp)) + { case AKAT: + tlo = des->cf[0]-pen(sp)*t[5]; + tup = des->cf[0]+pen(sp)*t[5]; + if ((tlo>cup) | (tup<clo)) /* done */ + i = 2; + else + { h1 = des->h; + clo = MAX(clo,tlo); + cup = MIN(cup,tup); + } + break; + case ACP: + cp = adcri(des->llk,t[0],t[2],pen(sp)); + if (cp<mcp) { mcp = cp; h1 = des->h; } + else + { if (i>0) i = 2; } + break; + case AMDI: + cp = mmse(lfd,sp,dv,des); + if (cp<mcp) { mcp = cp; h1 = des->h; } + else + { if (i>0) i = 2; } + } + } + return(h1); +} + +int alocfit(lfd,sp,dv,des,cv) +lfdata *lfd; +smpar *sp; +deriv *dv; +design *des; +int cv; +{ int lf_status; + double h0; + + lf_status = ainitband(lfd,sp,dv,des); + if (lf_error) return(lf_status); + if (acri(sp) == AOK) return(lf_status); + + h0 = fixh(sp); + fixh(sp) = aband2(lfd,sp,dv,des,des->h); + fixh(sp) = aband3(lfd,sp,dv,des,fixh(sp)); + nbhd(lfd,des,0,1,sp); + lf_status = locfit(lfd,des,sp,0,0,cv); + fixh(sp) = h0; + + return(lf_status); +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + * + * Evaluate the locfit fitting functions. + * calcp(sp,d) + * calculates the number of fitting functions. + * makecfn(sp,des,dv,d) + * makes the coef.number vector. + * fitfun(lfd, sp, x,t,f,dv) + * lfd is the local fit structure. + * sp smoothing parameter structure. + * x is the data point. + * t is the fitting point. + * f is a vector to return the results. + * dv derivative structure. + * designmatrix(lfd, sp, des) + * is a wrapper for fitfun to build the design matrix. + * + */ + +#include "locf.h" + +int calcp(sp,d) +smpar *sp; +int d; +{ int i, k; + + if (ubas(sp)) return(npar(sp)); + + switch (kt(sp)) + { case KSPH: + case KCE: + k = 1; + for (i=1; i<=deg(sp); i++) k = k*(d+i)/i; + return(k); + case KPROD: return(d*deg(sp)+1); + case KLM: return(d); + case KZEON: return(1); + } + LERR(("calcp: invalid kt %d",kt(sp))); + return(0); +} + +int coefnumber(dv,kt,d,deg) +int kt, d, deg; +deriv *dv; +{ int d0, d1, t; + + if (d==1) + { if (dv->nd<=deg) return(dv->nd); + return(-1); + } + + if (dv->nd==0) return(0); + if (deg==0) return(-1); + if (dv->nd==1) return(1+dv->deriv[0]); + if (deg==1) return(-1); + if (kt==KPROD) return(-1); + + if (dv->nd==2) + { d0 = dv->deriv[0]; d1 = dv->deriv[1]; + if (d0<d1) { t = d0; d0 = d1; d1 = t; } + return((d+1)*(d0+1)-d0*(d0+3)/2+d1); + } + if (deg==2) return(-1); + + LERR(("coefnumber not programmed for nd>=3")); + return(-1); +} + +void makecfn(sp,des,dv,d) +smpar *sp; +design *des; +deriv *dv; +int d; +{ int i, nd; + + nd = dv->nd; + + des->cfn[0] = coefnumber(dv,kt(sp),d,deg(sp)); + des->ncoef = 1; + if (nd >= deg(sp)) return; + if (kt(sp)==KZEON) return; + + if (d>1) + { if (nd>=2) return; + if ((nd>=1) && (kt(sp)==KPROD)) return; + } + + dv->nd = nd+1; + for (i=0; i<d; i++) + { dv->deriv[nd] = i; + des->cfn[i+1] = coefnumber(dv,kt(sp),d,deg(sp)); + } + dv->nd = nd; + + des->ncoef = 1+d; +} + +void fitfunangl(dx,ff,sca,cd,deg) +double dx, *ff, sca; +int deg, cd; +{ + if (deg>=3) WARN(("Can't handle angular model with deg>=3")); + + switch(cd) + { case 0: + ff[0] = 1; + ff[1] = sin(dx/sca)*sca; + ff[2] = (1-cos(dx/sca))*sca*sca; + return; + case 1: + ff[0] = 0; + ff[1] = cos(dx/sca); + ff[2] = sin(dx/sca)*sca; + return; + case 2: + ff[0] = 0; + ff[1] = -sin(dx/sca)/sca; + ff[2] = cos(dx/sca); + return; + default: WARN(("Can't handle angular model with >2 derivs")); + } +} + +void fitfun(lfd,sp,x,t,f,dv) +lfdata *lfd; +smpar *sp; +double *x, *t, *f; +deriv *dv; +{ int d, deg, nd, m, i, j, k, ct_deriv[MXDIM]; + double ff[MXDIM][1+MXDEG], dx[MXDIM], *xx[MXDIM]; + + if (ubas(sp)) + { for (i=0; i<lfd->d; i++) xx[i] = &x[i]; + i = 0; + sp->vbasis(xx,t,1,lfd->d,1,npar(sp),f); + return; + } + + d = lfd->d; + deg = deg(sp); + m = 0; + nd = (dv==NULL) ? 0 : dv->nd; + + if (kt(sp)==KZEON) + { f[0] = 1.0; + return; + } + + if (kt(sp)==KLM) + { for (i=0; i<d; i++) f[m++] = x[i]; + return; + } + + f[m++] = (nd==0); + if (deg==0) return; + + for (i=0; i<d; i++) + { ct_deriv[i] = 0; + dx[i] = (t==NULL) ? x[i] : x[i]-t[i]; + } + for (i=0; i<nd; i++) ct_deriv[dv->deriv[i]]++; + + for (i=0; i<d; i++) + { switch(lfd->sty[i]) + { + case STANGL: + fitfunangl(dx[i],ff[i],lfd->sca[i],ct_deriv[i],deg(sp)); + break; + default: + for (j=0; j<ct_deriv[i]; j++) ff[i][j] = 0.0; + ff[i][ct_deriv[i]] = 1.0; + for (j=ct_deriv[i]+1; j<=deg; j++) + ff[i][j] = ff[i][j-1]*dx[i]/(j-ct_deriv[i]); + } + } + +/* + * Product kernels. Note that if ct_deriv[i] != nd, that implies + * there is differentiation wrt another variable, and all components + * involving x[i] are 0. + */ + if ((d==1) || (kt(sp)==KPROD)) + { for (j=1; j<=deg; j++) + for (i=0; i<d; i++) + f[m++] = (ct_deriv[i]==nd) ? ff[i][j] : 0.0; + return; + } + +/* + * Spherical kernels with the full polynomial basis. + * Presently implemented up to deg=3. + */ + for (i=0; i<d; i++) + f[m++] = (ct_deriv[i]==nd) ? ff[i][1] : 0.0; + if (deg==1) return; + + for (i=0; i<d; i++) + { + /* xi^2/2 terms. */ + f[m++] = (ct_deriv[i]==nd) ? ff[i][2] : 0.0; + + /* xi xj terms */ + for (j=i+1; j<d; j++) + f[m++] = (ct_deriv[i]+ct_deriv[j]==nd) ? ff[i][1]*ff[j][1] : 0.0; + } + if (deg==2) return; + + for (i=0; i<d; i++) + { + /* xi^3/6 terms */ + f[m++] = (ct_deriv[i]==nd) ? ff[i][3] : 0.0; + + /* xi^2/2 xk terms */ + for (k=i+1; k<d; k++) + f[m++] = (ct_deriv[i]+ct_deriv[k]==nd) ? ff[i][2]*ff[k][1] : 0.0; + + /* xi xj xk terms */ + for (j=i+1; j<d; j++) + { f[m++] = (ct_deriv[i]+ct_deriv[j]==nd) ? ff[i][1]*ff[j][2] : 0.0; + for (k=j+1; k<d; k++) + f[m++] = (ct_deriv[i]+ct_deriv[j]+ct_deriv[k]==nd) ? + ff[i][1]*ff[j][1]*ff[k][1] : 0.0; + } + } + if (deg==3) return; + + LERR(("fitfun: can't handle deg=%d for spherical kernels",deg)); +} + +/* + * Build the design matrix. Assumes des->ind contains the indices of + * the required data points; des->n the number of points; des->xev + * the fitting point. + */ +void designmatrix(lfd,sp,des) +lfdata *lfd; +smpar *sp; +design *des; +{ int i, ii, j, p; + double *X, u[MXDIM]; + + X = d_x(des); + p = des->p; + + if (ubas(sp)) + { + sp->vbasis(lfd->x,des->xev,lfd->n,lfd->d,des->n,p,X); + return; + } + + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + for (j=0; j<lfd->d; j++) u[j] = datum(lfd,j,ii); + fitfun(lfd,sp,u,des->xev,&X[ii*p],NULL); + } +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + * + * + * Functions for determining bandwidth; smoothing neighborhood + * and smoothing weights. + */ + +#include "locf.h" + +double rho(x,sc,d,kt,sty) /* ||x|| for appropriate distance metric */ +double *x, *sc; +int d, kt, *sty; +{ double rhoi[MXDIM], s; + int i; + for (i=0; i<d; i++) + { if (sty!=NULL) + { switch(sty[i]) + { case STANGL: rhoi[i] = 2*sin(x[i]/(2*sc[i])); break; + case STCPAR: rhoi[i] = 0; break; + default: rhoi[i] = x[i]/sc[i]; + } } + else rhoi[i] = x[i]/sc[i]; + } + + if (d==1) return(fabs(rhoi[0])); + + s = 0; + if (kt==KPROD) + { for (i=0; i<d; i++) + { rhoi[i] = fabs(rhoi[i]); + if (rhoi[i]>s) s = rhoi[i]; + } + return(s); + } + + if (kt==KSPH) + { for (i=0; i<d; i++) + s += rhoi[i]*rhoi[i]; + return(sqrt(s)); + } + + LERR(("rho: invalid kt")); + return(0.0); +} + +double kordstat(x,k,n,ind) +double *x; +int k, n, *ind; +{ int i, i0, i1, l, r; + double piv; + if (k<1) return(0.0); + i0 = 0; i1 = n-1; + while (1) + { piv = x[ind[(i0+i1)/2]]; + l = i0; r = i1; + while (l<=r) + { while ((l<=i1) && (x[ind[l]]<=piv)) l++; + while ((r>=i0) && (x[ind[r]]>piv)) r--; + if (l<=r) ISWAP(ind[l],ind[r]); + } /* now, x[ind[i0..r]] <= piv < x[ind[l..i1]] */ + if (r<k-1) i0 = l; /* go right */ + else /* put pivots in middle */ + { for (i=i0; i<=r; ) + if (x[ind[i]]==piv) { ISWAP(ind[i],ind[r]); r--; } + else i++; + if (r<k-1) return(piv); + i1 = r; + } + } +} + +/* check if i'th data point is in limits */ +int inlim(lfd,i) +lfdata *lfd; +int i; +{ int d, j, k; + double *xlim; + + xlim = lfd->xl; + d = lfd->d; + k = 1; + for (j=0; j<d; j++) + { if (xlim[j]<xlim[j+d]) + k &= ((datum(lfd,j,i)>=xlim[j]) & (datum(lfd,j,i)<=xlim[j+d])); + } + return(k); +} + +double compbandwid(di,ind,x,n,d,nn,fxh) +double *di, *x, fxh; +int n, d, nn, *ind; +{ int i; + double nnh; + + if (nn==0) return(fxh); + + if (nn<n) + nnh = kordstat(di,nn,n,ind); + else + { nnh = 0; + for (i=0; i<n; i++) nnh = MAX(nnh,di[i]); + nnh = nnh*exp(log(1.0*nn/n)/d); + } + return(MAX(fxh,nnh)); +} + +/* + fast version of nbhd for ordered 1-d data +*/ +void nbhd1(lfd,sp,des,k) +lfdata *lfd; +smpar *sp; +design *des; +int k; +{ double x, h, *xd, sc; + int i, l, r, m, n, z; + + n = lfd->n; + x = des->xev[0]; + xd = dvari(lfd,0); + sc = lfd->sca[0]; + + /* find closest data point to x */ + if (x<=xd[0]) z = 0; + else + if (x>=xd[n-1]) z = n-1; + else + { l = 0; r = n-1; + while (r-l>1) + { z = (r+l)/2; + if (xd[z]>x) r = z; + else l = z; + } + /* now, xd[0..l] <= x < x[r..n-1] */ + if ((x-xd[l])>(xd[r]-x)) z = r; else z = l; + } + /* closest point to x is xd[z] */ + + if (nn(sp)<0) /* user bandwidth */ + h = sp->vb(des->xev); + else + { if (k>0) /* set h to nearest neighbor bandwidth */ + { l = r = z; + if (l==0) r = k-1; + if (r==n-1) l = n-k; + while (r-l<k-1) + { if ((x-xd[l-1])<(xd[r+1]-x)) l--; else r++; + if (l==0) r = k-1; + if (r==n-1) l = n-k; + } + h = x-xd[l]; + if (h<xd[r]-x) h = xd[r]-x; + } + else h = 0; + h /= sc; + if (h<fixh(sp)) h = fixh(sp); + } + + m = 0; + if (xd[z]>x) z--; /* so xd[z]<=x */ + /* look left */ + for (i=z; i>=0; i--) if (inlim(lfd,i)) + { dist(des,i) = (x-xd[i])/sc; + wght(des,i) = weight(lfd, sp, &xd[i], &x, h, 1, dist(des,i)); + if (wght(des,i)>0) + { des->ind[m] = i; + m++; + } else i = 0; + } + /* look right */ + for (i=z+1; i<n; i++) if (inlim(lfd,i)) + { dist(des,i) = (xd[i]-x)/sc; + wght(des,i) = weight(lfd, sp, &xd[i], &x, h, 1, dist(des,i)); + if (wght(des,i)>0) + { des->ind[m] = i; + m++; + } else i = n; + } + + des->n = m; + des->h = h; +} + +void nbhd_zeon(lfd,des) +lfdata *lfd; +design *des; +{ int i, j, m, eq; + + m = 0; + for (i=0; i<lfd->n; i++) + { eq = 1; + for (j=0; j<lfd->d; j++) eq = eq && (des->xev[j] == datum(lfd,j,i)); + if (eq) + { wght(des,i) = 1; + des->ind[m] = i; + m++; + } + } + des->n = m; + des->h = 1.0; +} + +void nbhd(lfd,des,nn,redo,sp) +lfdata *lfd; +design *des; +int redo, nn; +smpar *sp; +{ int d, i, j, m, n; + double h, u[MXDIM]; + + if (lf_debug>1) mut_printf("nbhd: nn %d fixh %8.5f\n",nn,fixh(sp)); + + d = lfd->d; n = lfd->n; + + if (ker(sp)==WPARM) + { for (i=0; i<n; i++) + { wght(des,i) = 1.0; + des->ind[i] = i; + } + des->n = n; + return; + } + + if (kt(sp)==KZEON) + { nbhd_zeon(lfd,des); + return; + } + + if (kt(sp)==KCE) + { des->h = 0.0; + return; + } + + /* ordered 1-dim; use fast searches */ + if ((nn<=n) & (lfd->ord) & (ker(sp)!=WMINM) & (lfd->sty[0]!=STANGL)) + { nbhd1(lfd,sp,des,nn); + return; + } + + if (!redo) + { for (i=0; i<n; i++) + { for (j=0; j<d; j++) u[j] = datum(lfd,j,i)-des->xev[j]; + dist(des,i) = rho(u,lfd->sca,d,kt(sp),lfd->sty); + des->ind[i] = i; + } + } + else + for (i=0; i<n; i++) des->ind[i] = i; + + if (ker(sp)==WMINM) + { des->h = minmax(lfd,des,sp); + return; + } + + if (nn<0) + h = sp->vb(des->xev); + else + h = compbandwid(des->di,des->ind,des->xev,n,lfd->d,nn,fixh(sp)); + m = 0; + for (i=0; i<n; i++) if (inlim(lfd,i)) + { for (j=0; j<d; j++) u[j] = datum(lfd,j,i); + wght(des,i) = weight(lfd, sp, u, des->xev, h, 1, dist(des,i)); + if (wght(des,i)>0) + { des->ind[m] = i; + m++; + } + } + des->n = m; + des->h = h; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + * + * This file includes functions to solve for the scale estimate in + * local robust regression and likelihood. The main entry point is + * lf_robust(lfd,sp,des,mxit), + * called from the locfit() function. + * + * The update_rs(x) accepts a residual scale x as the argument (actually, + * it works on the log-scale). The function computes the local fit + * assuming this residual scale, and re-estimates the scale from this + * new fit. The final solution satisfies the fixed point equation + * update_rs(x)=x. The function lf_robust() automatically calls + * update_rs() through the fixed point iterations. + * + * The estimation of the scale from the fit is based on the sqrt of + * the median deviance of observations with non-zero weights (in the + * gaussian case, this is the median absolute residual). + * + * TODO: + * Should use smoothing weights in the median. + */ + +#include "locf.h" + +extern int lf_status; +double robscale; + +static lfdata *rob_lfd; +static smpar *rob_sp; +static design *rob_des; +static int rob_mxit; + +double median(x,n) +double *x; +int n; +{ int i, j, lt, eq, gt; + double lo, hi, s; + lo = hi = x[0]; + for (i=0; i<n; i++) + { lo = MIN(lo,x[i]); + hi = MAX(hi,x[i]); + } + if (lo==hi) return(lo); + lo -= (hi-lo); + hi += (hi-lo); + for (i=0; i<n; i++) + { if ((x[i]>lo) & (x[i]<hi)) + { s = x[i]; lt = eq = gt = 0; + for (j=0; j<n; j++) + { lt += (x[j]<s); + eq += (x[j]==s); + gt += (x[j]>s); + } + if ((2*(lt+eq)>n) && (2*(gt+eq)>n)) return(s); + if (2*(lt+eq)<=n) lo = s; + if (2*(gt+eq)<=n) hi = s; + } + } + return((hi+lo)/2); +} + +double nrobustscale(lfd,sp,des,rs) +lfdata *lfd; +smpar *sp; +design *des; +double rs; +{ int i, ii, p; + double link[LLEN], sc, sd, sw, e; + p = des->p; sc = sd = sw = 0.0; + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + fitv(des,ii) = base(lfd,ii)+innerprod(des->cf,d_xi(des,ii),p); + e = resp(lfd,ii)-fitv(des,ii); + stdlinks(link,lfd,sp,ii,fitv(des,ii),rs); + sc += wght(des,ii)*e*link[ZDLL]; + sd += wght(des,ii)*e*e*link[ZDDLL]; + sw += wght(des,ii); + } + + /* newton-raphson iteration for log(s) + -psi(ei/s) - log(s); s = e^{-th} + */ + rs *= exp((sc-sw)/(sd+sc)); + return(rs); +} + +double robustscale(lfd,sp,des) +lfdata *lfd; +smpar *sp; +design *des; +{ int i, ii, p, fam, lin, or; + double rs, link[LLEN]; + p = des->p; + fam = fam(sp); + lin = link(sp); + or = fami(sp)->robust; + fami(sp)->robust = 0; + + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + fitv(des,ii) = base(lfd,ii) + innerprod(des->cf,d_xi(des,ii),p); + links(fitv(des,ii),resp(lfd,ii),fami(sp),lin,link,cens(lfd,ii),prwt(lfd,ii),1.0); + des->res[i] = -2*link[ZLIK]; + } + fami(sp)->robust = or; + rs = sqrt(median(des->res,des->n)); + + if (rs==0.0) rs = 1.0; + return(rs); +} + +double update_rs(x) +double x; +{ double nx; + if (lf_status != LF_OK) return(x); + robscale = exp(x); + lfiter(rob_lfd,rob_sp,rob_des,rob_mxit); + if (lf_status != LF_OK) return(x); + + nx = log(robustscale(rob_lfd,rob_sp,rob_des)); + if (nx<x-0.2) nx = x-0.2; + return(nx); +} + +void lf_robust(lfd,sp,des,mxit) +lfdata *lfd; +design *des; +smpar *sp; +int mxit; +{ double x; + rob_lfd = lfd; + rob_des = des; + rob_sp = sp; + rob_mxit = mxit; + lf_status = LF_OK; + + x = log(robustscale(lfd,sp,des)); + + solve_fp(update_rs, x, 1.0e-6, mxit); +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + * Post-fitting functions to compute the local variance and + * influence functions. Also the local degrees of freedom + * calculations for adaptive smoothing. + */ + +#include "locf.h" + +extern double robscale; + +/* + vmat() computes (after the local fit..) the matrix + M2 = X^T W^2 V X. + M12 = (X^T W V X)^{-1} M2 + Also, for convenience, tr[0] = sum(wi) tr[1] = sum(wi^2). +*/ +void vmat(lfd, sp, des, M12, M2) +lfdata *lfd; +smpar *sp; +design *des; +double *M12, *M2; +{ int i, ii, p, nk, ok; + double link[LLEN], h, ww, tr0, tr1; + p = des->p; + setzero(M2,p*p); + + nk = -1; + + /* for density estimation, use integral rather than + sum form, if W^2 is programmed... + */ + if ((fam(sp)<=THAZ) && (link(sp)==LLOG)) + { switch(ker(sp)) + { case WGAUS: nk = WGAUS; h = des->h/SQRT2; break; + case WRECT: nk = WRECT; h = des->h; break; + case WEPAN: nk = WBISQ; h = des->h; break; + case WBISQ: nk = WQUQU; h = des->h; break; + case WTCUB: nk = W6CUB; h = des->h; break; + case WEXPL: nk = WEXPL; h = des->h/2; break; + } + } + + tr0 = tr1 = 0.0; + if (nk != -1) + { ok = ker(sp); ker(sp) = nk; +/* compute M2 using integration. Use M12 as work matrix. */ + (des->itype)(des->xev, M2, M12, des->cf, h); + ker(sp) = ok; + if (fam(sp)==TDEN) multmatscal(M2,des->smwt,p*p); + tr0 = des->ss[0]; + tr1 = M2[0]; /* n int W e^<a,A> */ + } + else + { for (i=0; i<des->n; i++) + { ii = des->ind[i]; + stdlinks(link,lfd,sp,ii,fitv(des,ii),robscale); + ww = SQR(wght(des,ii))*link[ZDDLL]; + tr0 += wght(des,ii); + tr1 += SQR(wght(des,ii)); + addouter(M2,d_xi(des,ii),d_xi(des,ii),p,ww); + } + } + des->tr0 = tr0; + des->tr1 = tr1; + + memcpy(M12,M2,p*p*sizeof(double)); + for (i=0; i<p; i++) + jacob_solve(&des->xtwx,&M12[i*p]); +} + +void lf_vcov(lfd,sp,des) +lfdata *lfd; +smpar *sp; +design *des; +{ int i, j, k, p; + double *M12, *M2; + M12 = des->V; M2 = des->P; p = des->p; + vmat(lfd,sp,des,M12,M2); /* M2 = X^T W^2 V X tr0=sum(W) tr1=sum(W*W) */ + des->tr2 = m_trace(M12,p); /* tr (XTWVX)^{-1}(XTW^2VX) */ + +/* + * Covariance matrix is M1^{-1} * M2 * M1^{-1} + * We compute this using the cholesky decomposition of + * M2; premultiplying by M1^{-1} and squaring. This + * is more stable than direct computation in near-singular cases. + */ + chol_dec(M2,p,p); + for (i=0; i<p; i++) + for (j=0; j<i; j++) + { M2[j*p+i] = M2[i*p+j]; + M2[i*p+j] = 0.0; + } + for (i=0; i<p; i++) jacob_solve(&des->xtwx,&M2[i*p]); + for (i=0; i<p; i++) + { for (j=0; j<p; j++) + { M12[i*p+j] = 0; + for (k=0; k<p; k++) + M12[i*p+j] += M2[k*p+i]*M2[k*p+j]; /* ith column of covariance */ + } + } + if ((fam(sp)==TDEN) && (link(sp)==LIDENT)) + multmatscal(M12,1/SQR(des->smwt),p*p); + +/* this computes the influence function as des->f1[0]. */ + unitvec(des->f1,0,des->p); + jacob_solve(&des->xtwx,des->f1); +} + +/* local_df computes: + * tr[0] = trace(W) + * tr[1] = trace(W*W) + * tr[2] = trace( M1^{-1} M2 ) + * tr[3] = trace( M1^{-1} M3 ) + * tr[4] = trace( (M1^{-1} M2)^2 ) + * tr[5] = var(theta-hat). + */ +void local_df(lfd,sp,des,tr) +lfdata *lfd; +smpar *sp; +design *des; +double *tr; +{ int i, ii, j, p; + double *m2, *V, ww, link[LLEN]; + + tr[0] = tr[1] = tr[2] = tr[3] = tr[4] = tr[5] = 0.0; + m2 = des->V; V = des->P; p = des->p; + + vmat(lfd,sp,des,m2,V); /* M = X^T W^2 V X tr0=sum(W) tr1=sum(W*W) */ + tr[0] = des->tr0; + tr[1] = des->tr1; + tr[2] = m_trace(m2,p); /* tr (XTWVX)^{-1}(XTW^2VX) */ + + unitvec(des->f1,0,p); + jacob_solve(&des->xtwx,des->f1); + for (i=0; i<p; i++) + for (j=0; j<p; j++) + { tr[4] += m2[i*p+j]*m2[j*p+i]; /* tr(M^2) */ + tr[5] += des->f1[i]*V[i*p+j]*des->f1[j]; /* var(thetahat) */ + } + tr[5] = sqrt(tr[5]); + + setzero(m2,p*p); + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + stdlinks(link,lfd,sp,ii,fitv(des,ii),robscale); + ww = wght(des,ii)*wght(des,ii)*wght(des,ii)*link[ZDDLL]; + addouter(m2,d_xi(des,ii),d_xi(des,ii),p,ww); + } + for (i=0; i<p; i++) + { jacob_solve(&des->xtwx,&m2[i*p]); + tr[3] += m2[i*(p+1)]; + } + + return; +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + * Routines for computing weight diagrams. + * wdiag(lf,des,lx,deg,ty,exp) + * Must locfit() first, unless ker==WPARM and has par. comp. + * + */ + +#include "locf.h" + +static double *wd; +extern double robscale; +void nnresproj(lfd,sp,des,u,m,p) +lfdata *lfd; +smpar *sp; +design *des; +double *u; +int m, p; +{ int i, ii, j; + double link[LLEN]; + setzero(des->f1,p); + for (j=0; j<m; j++) + { ii = des->ind[j]; + stdlinks(link,lfd,sp,ii,fitv(des,ii),robscale); + for (i=0; i<p; i++) des->f1[i] += link[ZDDLL]*d_xij(des,j,ii)*u[j]; + } + jacob_solve(&des->xtwx,des->f1); + for (i=0; i<m; i++) + { ii = des->ind[i]; + u[i] -= innerprod(des->f1,d_xi(des,ii),p)*wght(des,ii); + } +} + +void wdexpand(l,n,ind,m) +double *l; +int *ind, n, m; +{ int i, j, t; + double z; + for (j=m; j<n; j++) { l[j] = 0.0; ind[j] = -1; } + j = m-1; + while (j>=0) + { if (ind[j]==j) j--; + else + { i = ind[j]; + z = l[j]; l[j] = l[i]; l[i] = z; + t = ind[j]; ind[j] = ind[i]; ind[i] = t; + if (ind[j]==-1) j--; + } + } + +/* for (i=n-1; i>=0; i--) + { l[i] = ((j>=0) && (ind[j]==i)) ? l[j--] : 0.0; } */ +} + +int wdiagp(lfd,sp,des,lx,pc,dv,deg,ty,exp) +lfdata *lfd; +smpar *sp; +design *des; +paramcomp *pc; +deriv *dv; +double *lx; +int deg, ty, exp; +{ int i, j, p, nd; + double *l1; + + p = des->p; + + fitfun(lfd,sp,des->xev,pc->xbar,des->f1,dv); + if (exp) + { jacob_solve(&pc->xtwx,des->f1); + for (i=0; i<lfd->n; i++) + lx[i] = innerprod(des->f1,d_xi(des,des->ind[i]),p); + return(lfd->n); + } + jacob_hsolve(&pc->xtwx,des->f1); + for (i=0; i<p; i++) lx[i] = des->f1[i]; + + nd = dv->nd; + dv->nd = nd+1; + if (deg>=1) + for (i=0; i<lfd->d; i++) + { dv->deriv[nd] = i; + l1 = &lx[(i+1)*p]; + fitfun(lfd,sp,des->xev,pc->xbar,l1,dv); + jacob_hsolve(&pc->xtwx,l1); + } + + dv->nd = nd+2; + if (deg>=2) + for (i=0; i<lfd->d; i++) + { dv->deriv[nd] = i; + for (j=0; j<lfd->d; j++) + { dv->deriv[nd+1] = j; + l1 = &lx[(i*lfd->d+j+lfd->d+1)*p]; + fitfun(lfd,sp,des->xev,pc->xbar,l1,dv); + jacob_hsolve(&pc->xtwx,l1); + } } + dv->nd = nd; + return(p); +} + +int wdiag(lfd,sp,des,lx,dv,deg,ty,exp) +lfdata *lfd; +smpar *sp; +design *des; +deriv *dv; +double *lx; +int deg, ty, exp; +/* deg=0: l(x) only. + deg=1: l(x), l'(x) + deg=2: l(x), l'(x), l''(x) + ty = 1: e1 (X^T WVX)^{-1} X^T W -- hat matrix + ty = 2: e1 (X^T WVX)^{-1} X^T WV^{1/2} -- scb's +*/ +{ double w, *X, *lxd, *lxdd, wdd, wdw, *ulx, link[LLEN], h; + double dfx[MXDIM], hs[MXDIM]; + int i, ii, j, k, l, m, d, p, nd; + + h = des->h; + nd = dv->nd; + wd = des->wd; + d = lfd->d; p = des->p; X = d_x(des); + ulx = des->res; + m = des->n; + for (i=0; i<d; i++) hs[i] = h*lfd->sca[i]; + if (deg>0) + { lxd = &lx[m]; + setzero(lxd,m*d); + if (deg>1) + { lxdd = &lxd[d*m]; + setzero(lxdd,m*d*d); + } } + + if (nd>0) fitfun(lfd,sp,des->xev,des->xev,des->f1,dv); /* c(0) */ + else unitvec(des->f1,0,p); + jacob_solve(&des->xtwx,des->f1); /* c(0) (X^TWX)^{-1} */ + for (i=0; i<m; i++) + { ii = des->ind[i]; + lx[i] = innerprod(des->f1,&X[ii*p],p); /* c(0)(XTWX)^{-1}X^T */ + if (deg>0) + { wd[i] = Wd(dist(des,ii)/h,ker(sp)); + for (j=0; j<d; j++) + { dfx[j] = datum(lfd,j,ii)-des->xev[j]; + lxd[j*m+i] = lx[i]*wght(des,ii)*weightd(dfx[j],lfd->sca[j], + d,ker(sp),kt(sp),h,lfd->sty[j],dist(des,ii)); + /* c(0) (XTWX)^{-1}XTW' */ + } + if (deg>1) + { wdd = Wdd(dist(des,ii)/h,ker(sp)); + for (j=0; j<d; j++) + for (k=0; k<d; k++) + { w = (dist(des,ii)==0) ? 0 : h/dist(des,ii); + w = wdd * (des->xev[k]-datum(lfd,k,ii)) * (des->xev[j]-datum(lfd,j,ii)) + * w*w / (hs[k]*hs[k]*hs[j]*hs[j]); + if (j==k) w += wd[i]/(hs[j]*hs[j]); + lxdd[(j*d+k)*m+i] = lx[i]*w; + /* c(0)(XTWX)^{-1}XTW'' */ + } + } + } + lx[i] *= wght(des,ii); + } + + dv->nd = nd+1; + if (deg==2) + { for (i=0; i<d; i++) + { dv->deriv[nd] = i; + fitfun(lfd,sp,des->xev,des->xev,des->f1,dv); + for (k=0; k<m; k++) + { ii = des->ind[i]; + stdlinks(link,lfd,sp,ii,fitv(des,ii),robscale); + for (j=0; j<p; j++) + des->f1[j] -= link[ZDDLL]*lxd[i*m+k]*X[ii*p+j]; + /* c'(x)-c(x)(XTWX)^{-1}XTW'X */ + } + jacob_solve(&des->xtwx,des->f1); /* (...)(XTWX)^{-1} */ + for (j=0; j<m; j++) + { ii = des->ind[j]; + ulx[j] = innerprod(des->f1,&X[ii*p],p); /* (...)XT */ + } + for (j=0; j<d; j++) + for (k=0; k<m; k++) + { ii = des->ind[k]; + dfx[j] = datum(lfd,j,ii)-des->xev[j]; + wdw = wght(des,ii)*weightd(dfx[j],lfd->sca[j],d,ker(sp), + kt(sp),h,lfd->sty[j],dist(des,ii)); + lxdd[(i*d+j)*m+k] += ulx[k]*wdw; + lxdd[(j*d+i)*m+k] += ulx[k]*wdw; + } /* + 2(c'-c(XTWX)^{-1}XTW'X)(XTWX)^{-1}XTW' */ + } + for (j=0; j<d*d; j++) nnresproj(lfd,sp,des,&lxdd[j*m],m,p); + /* * (I-X(XTWX)^{-1} XTW */ + } + if (deg>0) + { for (j=0; j<d; j++) nnresproj(lfd,sp,des,&lxd[j*m],m,p); + /* c(0)(XTWX)^{-1}XTW'(I-X(XTWX)^{-1}XTW) */ + for (i=0; i<d; i++) + { dv->deriv[nd]=i; + fitfun(lfd,sp,des->xev,des->xev,des->f1,dv); + jacob_solve(&des->xtwx,des->f1); + for (k=0; k<m; k++) + { ii = des->ind[k]; + for (l=0; l<p; l++) + lxd[i*m+k] += des->f1[l]*X[ii*p+l]*wght(des,ii); + } /* add c'(0)(XTWX)^{-1}XTW */ + } + } + + dv->nd = nd+2; + if (deg==2) + { for (i=0; i<d; i++) + { dv->deriv[nd]=i; + for (j=0; j<d; j++) + { dv->deriv[nd+1]=j; + fitfun(lfd,sp,des->xev,des->xev,des->f1,dv); + jacob_solve(&des->xtwx,des->f1); + for (k=0; k<m; k++) + { ii = des->ind[k]; + for (l=0; l<p; l++) + lxdd[(i*d+j)*m+k] += des->f1[l]*X[ii*p+l]*wght(des,ii); + } /* + c''(x)(XTWX)^{-1}XTW */ + } + } + } + dv->nd = nd; + + k = 1+d*(deg>0)+d*d*(deg==2); + + if (exp) wdexpand(lx,lfd->n,des->ind,m); + + if (ty==1) return(m); + for (i=0; i<m; i++) + { ii = des->ind[i]; + stdlinks(link,lfd,sp,ii,fitv(des,ii),robscale); + link[ZDDLL] = sqrt(fabs(link[ZDDLL])); + for (j=0; j<k; j++) lx[j*m+i] *= link[ZDDLL]; + } + return(m); +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + * String matching functions. For a given argument string, find + * the best match from an array of possibilities. Is there a library + * function somewhere to do something like this? + * + * return values of -1 indicate failure/unknown string. + */ + +#include "locf.h" + +int ct_match(z1, z2) +char *z1, *z2; +{ int ct = 0; + while (z1[ct]==z2[ct]) + { if (z1[ct]=='\0') return(ct+1); + ct++; + } + return(ct); +} + +int pmatch(z, strings, vals, n, def) +char *z, **strings; +int *vals, n, def; +{ int i, ct, best, best_ct; + best = -1; + best_ct = 0; + + for (i=0; i<n; i++) + { ct = ct_match(z,strings[i]); + if (ct==strlen(z)+1) return(vals[i]); + if (ct>best_ct) { best = i; best_ct = ct; } + } + if (best==-1) return(def); + return(vals[best]); +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +#include "locf.h" + +int lf_maxit = 20; +int lf_debug = 0; +int lf_error = 0; + +double s0, s1; +static lfdata *lf_lfd; +static design *lf_des; +static smpar *lf_sp; +int lf_status; +int ident=0; +double lf_tol; +extern double robscale; + +void lfdata_init(lfd) +lfdata *lfd; +{ int i; + for (i=0; i<MXDIM; i++) + { lfd->sty[i] = 0; + lfd->sca[i] = 1.0; + lfd->xl[i] = lfd->xl[i+MXDIM] = 0.0; + } + lfd->y = lfd->w = lfd->c = lfd->b = NULL; + lfd->d = lfd->n = 0; +} + +void smpar_init(sp,lfd) +smpar *sp; +lfdata *lfd; +{ nn(sp) = 0.7; + fixh(sp)= 0.0; + pen(sp) = 0.0; + acri(sp)= ANONE; + deg(sp) = deg0(sp) = 2; + ubas(sp) = 0; + kt(sp) = KSPH; + ker(sp) = WTCUB; + fam(sp) = 64+TGAUS; + link(sp)= LDEFAU; + npar(sp) = calcp(sp,lfd->d); +} + +void deriv_init(dv) +deriv *dv; +{ dv->nd = 0; +} + +int des_reqd(n,p) +int n, p; +{ + return(n*(p+5)+2*p*p+4*p + jac_reqd(p)); +} +int des_reqi(n,p) +int n, p; +{ return(n+p); +} + +void des_init(des,n,p) +design *des; +int n, p; +{ double *z; + int k; + + if (n<=0) WARN(("des_init: n <= 0")); + if (p<=0) WARN(("des_init: p <= 0")); + + if (des->des_init_id != DES_INIT_ID) + { des->lwk = des->lind = 0; + des->des_init_id = DES_INIT_ID; + } + + k = des_reqd(n,p); + if (k>des->lwk) + { des->wk = (double *)calloc(k,sizeof(double)); + if ( des->wk == NULL ) { + printf("Problem allocating memory for des->wk\n");fflush(stdout); + } + des->lwk = k; + } + z = des->wk; + + des->X = z; z += n*p; + des->w = z; z += n; + des->res=z; z += n; + des->di =z; z += n; + des->th =z; z += n; + des->wd =z; z += n; + des->V =z; z += p*p; + des->P =z; z += p*p; + des->f1 =z; z += p; + des->ss =z; z += p; + des->oc =z; z += p; + des->cf =z; z += p; + + z = jac_alloc(&des->xtwx,p,z); + + k = des_reqi(n,p); + if (k>des->lind) + { + des->ind = (int *)calloc(k,sizeof(int)); + if ( des->ind == NULL ) { + printf("Problem allocating memory for des->ind\n");fflush(stdout); + } + des->lind = k; + } + des->fix = &des->ind[n]; + for (k=0; k<p; k++) des->fix[k] = 0; + + des->n = n; des->p = p; + des->smwt = n; + des->xtwx.p = p; +} + +void deschk(des,n,p) +design *des; +int n, p; +{ WARN(("deschk deprecated - use des_init()")); + des_init(des,n,p); +} + +int likereg(coef, lk0, f1, Z) +double *coef, *lk0, *f1, *Z; +{ int i, ii, j, p; + double lk, ww, link[LLEN], *X; + + if (lf_debug>2) mut_printf(" likereg: %8.5f\n",coef[0]); + lf_status = LF_OK; + lk = 0.0; p = lf_des->p; + setzero(Z,p*p); + setzero(f1,p); + for (i=0; i<lf_des->n; i++) + { + ii = lf_des->ind[i]; + X = d_xi(lf_des,ii); + fitv(lf_des,ii) = base(lf_lfd,ii)+innerprod(coef,X,p); + lf_status = stdlinks(link,lf_lfd,lf_sp,ii,fitv(lf_des,ii),robscale); + if (lf_status == LF_BADP) + { *lk0 = -1.0e300; + return(NR_REDUCE); + } + if (lf_error) lf_status = LF_ERR; + if (lf_status != LF_OK) return(NR_BREAK); + + ww = wght(lf_des,ii); + lk += ww*link[ZLIK]; + for (j=0; j<p; j++) + f1[j] += X[j]*ww*link[ZDLL]; + addouter(Z, X, X, p, ww*link[ZDDLL]); + } + for (i=0; i<p; i++) if (lf_des->fix[i]) + { for (j=0; j<p; j++) Z[i*p+j] = Z[j*p+i] = 0.0; + Z[i*p+i] = 1.0; + f1[i] = 0.0; + } + + if (lf_debug>4) prresp(coef,Z,p); + if (lf_debug>3) mut_printf(" likelihood: %8.5f\n",lk); + *lk0 = lf_des->llk = lk; + + lf_status = fami(lf_sp)->pcheck(lf_sp,lf_des,lf_lfd); + switch(lf_status) + { case LF_DONE: return(NR_BREAK); + case LF_OOB: return(NR_REDUCE); + case LF_PF: return(NR_REDUCE); + case LF_NSLN: return(NR_BREAK); + } + + return(NR_OK); +} + +int reginit(lfd,des,sp) +lfdata *lfd; +design *des; +smpar *sp; +{ int i, ii; + double sb, link[LLEN]; + s0 = s1 = sb = 0; + for (i=0; i<des->n; i++) + { ii = des->ind[i]; + links(base(lfd,ii),resp(lfd,ii),fami(sp),LINIT,link,cens(lfd,ii),prwt(lfd,ii),1.0); + s1 += wght(des,ii)*link[ZDLL]; + s0 += wght(des,ii)*prwt(lfd,ii); + sb += wght(des,ii)*prwt(lfd,ii)*base(lfd,ii); + } + if (s0==0) return(LF_NOPT); /* no observations with W>0 */ + setzero(des->cf,des->p); + lf_tol = 1.0e-6*s0; + switch(link(sp)) + { case LIDENT: + des->cf[0] = (s1-sb)/s0; + return(LF_OK); + case LLOG: + if (s1<=0.0) + { des->cf[0] = -1000; + return(LF_INFA); + } + des->cf[0] = log(s1/s0) - sb/s0; + return(LF_OK); + case LLOGIT: + if (s1<=0.0) + { des->cf[0] = -1000; + return(LF_INFA); + } + if (s1>=s0) + { des->cf[0] = 1000; + return(LF_INFA); + } + des->cf[0] = logit(s1/s0)-sb/s0; + return(LF_OK); + case LINVER: + if (s1<=0.0) + { des->cf[0] = 1e100; + return(LF_INFA); + } + des->cf[0] = s0/s1-sb/s0; + return(LF_OK); + case LSQRT: + des->cf[0] = sqrt(s1/s0)-sb/s0; + return(LF_OK); + case LASIN: + des->cf[0] = asin(sqrt(s1/s0))-sb/s0; + return(LF_OK); + default: + LERR(("reginit: invalid link %d",link(sp))); + return(LF_ERR); + } +} + +int lfinit(lfd,sp,des) +lfdata *lfd; +smpar *sp; +design *des; +{ int initstat; + des->xtwx.sm = (deg0(sp)<deg(sp)) ? JAC_CHOL : JAC_EIGD; + + designmatrix(lfd,sp,des); + setfamily(sp); + initstat = fami(sp)->initial(lfd,des,sp); + + return(initstat); +} + +void lfiter(lfd,sp,des,maxit) +lfdata *lfd; +smpar *sp; +design *des; +int maxit; +{ int err; + if (lf_debug>1) mut_printf(" lfiter: %8.5f\n",des->cf[0]); + + lf_des = des; + lf_lfd = lfd; + lf_sp = sp; + + max_nr(fami(sp)->like, des->cf, des->oc, des->res, des->f1, + &des->xtwx, des->p, maxit, lf_tol, &err); + switch(err) + { case NR_OK: return; + case NR_NCON: + WARN(("max_nr not converged")); + return; + case NR_NDIV: + WARN(("max_nr reduction problem")); + return; + } + WARN(("max_nr return status %d",err)); +} + +int use_robust_scale(int tg) +{ if ((tg&64)==0) return(0); /* not quasi - no scale */ + if (((tg&128)==0) & (((tg&63)!=TROBT) & ((tg&63)!=TCAUC))) return(0); + return(1); +} + +/* + * noit not really needed any more, since + * gauss->pcheck returns LF_DONE, and likereg NR_BREAK + * in gaussian case. + * nb: 0/1: does local neighborhood and weights need computing? + * cv: 0/1: is variance/covariance matrix needed? + */ +int locfit(lfd,des,sp,noit,nb,cv) +lfdata *lfd; +design *des; +smpar *sp; +int noit, nb, cv; +{ int i; + + if (des->xev==NULL) + { LERR(("locfit: NULL evaluation point?")); + return(246); + } + + if (lf_debug>0) + { mut_printf("locfit: "); + for (i=0; i<lfd->d; i++) mut_printf(" %10.6f",des->xev[i]); + mut_printf("\n"); + } + +/* the 1e-12 avoids problems that can occur with roundoff */ + if (nb) nbhd(lfd,des,(int)(lfd->n*nn(sp)+1e-12),0,sp); + + lf_status = lfinit(lfd,sp,des); + + if (lf_status == LF_OK) + { if (use_robust_scale(fam(sp))) + lf_robust(lfd,sp,des,lf_maxit); + else + { if ((fam(sp)&63)==TQUANT) + lfquantile(lfd,sp,des,lf_maxit); + else + { robscale = 1.0; + lfiter(lfd,sp,des,lf_maxit); + } + } + } + + if (lf_status == LF_DONE) lf_status = LF_OK; + if (lf_status == LF_OOB) lf_status = LF_OK; + + if ((fam(sp)&63)==TDEN) /* convert from rate to density */ + { switch(link(sp)) + { case LLOG: + des->cf[0] -= log(des->smwt); + break; + case LIDENT: + multmatscal(des->cf,1.0/des->smwt,des->p); + break; + default: LERR(("Density adjustment; invalid link")); + } + } + + /* variance calculations, if requested */ + if (cv) + { switch(lf_status) + { case LF_PF: /* for these cases, variance calc. would likely fail. */ + case LF_NOPT: + case LF_NSLN: + case LF_INFA: + case LF_DEMP: + case LF_XOOR: + case LF_DNOP: + case LF_BADP: + des->llk = des->tr0 = des->tr1 = des->tr2 = 0.0; + setzero(des->V,des->p*des->p); + setzero(des->f1,des->p); + break; + default: lf_vcov(lfd,sp,des); + } + } + + return(lf_status); +} + +void lf_status_msg(status) +int status; +{ switch(status) +{ case LF_OK: return; + case LF_NCON: WARN(("locfit did not converge")); return; + case LF_OOB: WARN(("parameters out of bounds")); return; + case LF_PF: WARN(("perfect fit")); return; + case LF_NOPT: WARN(("no points with non-zero weight")); return; + case LF_NSLN: WARN(("no solution")); return; + case LF_INFA: WARN(("initial value problem")); return; + case LF_DEMP: WARN(("density estimate, empty integration region")); return; + case LF_XOOR: WARN(("procv: fit point outside xlim region")); return; + case LF_DNOP: WARN(("density estimation -- insufficient points in smoothing window")); return; + case LF_BADP: WARN(("bad parameters")); return; + default: WARN(("procv: unknown return code %d",status)); return; +} } +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + * Compute minimax weights for local regression. + */ + +#include "locf.h" +#define NR_EMPTY 834 + +int mmsm_ct; + +static int debug=0; +#define CONVTOL 1.0e-8 +#define SINGTOL 1.0e-10 +#define NR_SINGULAR 100 + +static lfdata *mm_lfd; +static design *mm_des; +static double mm_gam, mmf, lb; +static int st; + +double ipower(x,n) /* use for n not too large!! */ +double x; +int n; +{ if (n==0) return(1.0); + if (n<0) return(1/ipower(x,-n)); + return(x*ipower(x,n-1)); +} + +double setmmwt(des,a,gam) +design *des; +double *a, gam; +{ double ip, w0, w1, sw, wt; + int i; + sw = 0.0; + for (i=0; i<mm_lfd->n; i++) + { ip = innerprod(a,d_xi(des,i),des->p); + wt = prwt(mm_lfd,i); + w0 = ip - gam*des->wd[i]; + w1 = ip + gam*des->wd[i]; + wght(des,i) = 0.0; + if (w0>0) { wght(des,i) = w0; sw += wt*w0*w0; } + if (w1<0) { wght(des,i) = w1; sw += wt*w1*w1; } + } + return(sw/2-a[0]); +} + +/* compute sum_{w!=0} AA^T; e1-sum wA */ +int mmsums(des,coef,f,z,J) +design *des; +double *coef, *f, *z; +jacobian *J; +{ int ct, i, j, p, sing; + double *A; + +mmsm_ct++; + A = J->Z; + *f = setmmwt(des,coef,mm_gam); + + p = des->p; + setzero(A,p*p); + setzero(z,p); + z[0] = 1.0; + ct = 0; + + for (i=0; i<mm_lfd->n; i++) + if (wght(des,i)!=0.0) + { addouter(A,d_xi(des,i),d_xi(des,i),p,prwt(mm_lfd,i)); + for (j=0; j<p; j++) z[j] -= prwt(mm_lfd,i)*wght(des,i)*d_xij(des,i,j); + ct++; + } + if (ct==0) return(NR_EMPTY); + + J->st = JAC_RAW; + J->p = p; + jacob_dec(J,JAC_EIGD); + + sing = 0; + for (i=0; i<p; i++) sing |= (J->Z[i*p+i]<SINGTOL); + if ((debug) & (sing)) mut_printf("SINGULAR!!!!\n"); + + return((sing) ? NR_SINGULAR : NR_OK); +} + +int descenddir(des,coef,dlt,f,af) +design *des; +double *coef, *dlt, *f; +int af; +{ int i, p; + double f0, *oc; + + if (debug) mut_printf("descenddir: %8.5f %8.5f\n",dlt[0],dlt[1]); + + f0 = *f; + oc = des->oc; + p = des->p; + memcpy(oc,coef,p*sizeof(double)); + + for (i=0; i<p; i++) coef[i] = oc[i]+lb*dlt[i]; + st = mmsums(des,coef,f,des->f1,&des->xtwx); + + if (*f>f0) /* halve till we drop */ + { while (*f>f0) + { lb = lb/2.0; + for (i=0; i<p; i++) coef[i] = oc[i]+lb*dlt[i]; + st = mmsums(des,coef,f,des->f1,&des->xtwx); + } + return(st); + } + + if (!af) return(st); + + /* double */ + while (*f<f0) + { f0 = *f; + lb *= 2.0; + for (i=0; i<p; i++) coef[i] = oc[i]+lb*dlt[i]; + st = mmsums(des,coef,f,des->f1,&des->xtwx); + } + + lb /= 2.0; + for (i=0; i<p; i++) coef[i] = oc[i]+lb*dlt[i]; + st = mmsums(des,coef,f,des->f1,&des->xtwx); + + return(st); +} + +int mm_initial(des) +design *des; +{ double *dlt; + + dlt = des->ss; + + setzero(des->cf,des->p); + st = mmsums(des,des->cf,&mmf,des->f1,&des->xtwx); + + setzero(dlt,des->p); + dlt[0] = 1; + lb = 1.0; + st = descenddir(des,des->cf,dlt,&mmf,1); + return(st); +} + +void getsingdir(des,dlt) +design *des; +double *dlt; +{ double f, sw, c0; + int i, j, p, sd; + + sd = -1; p = des->p; + setzero(dlt,p); + for (i=0; i<p; i++) if (des->xtwx.Z[i*p+i]<SINGTOL) sd = i; + if (sd==-1) + { mut_printf("getsingdir: nonsing?\n"); + return; + } + if (des->xtwx.dg[sd]>0) + for (i=0; i<p; i++) dlt[i] = des->xtwx.Q[p*i+sd]*des->xtwx.dg[i]; + else + { dlt[sd] = 1.0; + } + + c0 = innerprod(dlt,des->f1,p); + if (c0<0) for (i=0; i<p; i++) dlt[i] = -dlt[i]; +} + +void mmax(coef, old_coef, delta, J, p, maxit, tol, err) +double *coef, *old_coef, *delta, tol; +int p, maxit, *err; +jacobian *J; +{ double old_f, lambda; + int i, j; + + *err = NR_OK; + + for (j=0; j<maxit; j++) + { memcpy(old_coef,coef,p*sizeof(double)); + old_f = mmf; + + if (st == NR_SINGULAR) + { + getsingdir(mm_des,delta); + st = descenddir(mm_des,coef,delta,&mmf,1); + } + if (st == NR_EMPTY) + { + setzero(delta,p); + delta[0] = 1.0; + st = descenddir(mm_des,coef,delta,&mmf,1); + } + if (st == NR_OK) + { + lb = 1.0; + jacob_solve(J,mm_des->f1); + memcpy(delta,mm_des->f1,p*sizeof(double)); + st = descenddir(mm_des,coef,delta,&mmf,0); + } + + if ((j>0) & (fabs(mmf-old_f)<tol)) return; + } + WARN(("findab not converged")); + *err = NR_NCON; + return; +} + +double findab(gam) +double gam; +{ double sl; + int i, p, nr_stat; + + if (debug) mut_printf(" findab: gam %8.5f\n",gam); + mm_gam = gam; + p = mm_des->p; + lb = 1.0; + st = mm_initial(mm_des); + + mmax(mm_des->cf, mm_des->oc, mm_des->ss, + &mm_des->xtwx, p, lf_maxit, CONVTOL, &nr_stat); + + sl = 0.0; + for (i=0; i<mm_lfd->n; i++) sl += fabs(wght(mm_des,i))*mm_des->wd[i]; + + if (debug) mut_printf(" sl %8.5f gam %8.5f %8.5f %d\n", sl,gam,sl-gam,nr_stat); + return(sl-gam); +} + +double weightmm(coef,di,ff,gam) +double *coef, di, *ff, gam; +{ double y1, y2, ip; + ip = innerprod(ff,coef,mm_des->p); + y1 = ip-gam*di; if (y1>0) return(y1/ip); + y2 = ip+gam*di; if (y2<0) return(y2/ip); + return(0.0); +} + +double minmax(lfd,des,sp) +lfdata *lfd; +design *des; +smpar *sp; +{ double h, u[MXDIM], gam; + int i, j, m, d1, p1, err_flag; + + if (debug) mut_printf("minmax: x %8.5f\n",des->xev[0]); + mm_lfd = lfd; + mm_des = des; + +mmsm_ct = 0; + d1 = deg(sp)+1; + p1 = factorial(d1); + for (i=0; i<lfd->n; i++) + { for (j=0; j<lfd->d; j++) u[j] = datum(lfd,j,i); + des->wd[i] = sp->nn/p1*ipower(dist(des,i),d1); + des->ind[i] = i; + fitfun(lfd, sp, u,des->xev,d_xi(des,i),NULL); + } + +/* find gamma (i.e. solve eqn 13.17 from book), using the secant method. + * As a side effect, this finds the other minimax coefficients. + * Note that 13.17 is rewritten as + * g2 = sum |l_i(x)| (||xi-x||^(p+1) M/(s*(p+1)!)) + * where g2 = gamma * s * (p+1)! / M. The gam variable below is g2. + * The smoothing parameter is sp->nn == M/s. + */ + gam = solve_secant(findab, 0.0, 0.0,1.0, 0.0000001, BDF_EXPRIGHT, &err_flag); + +/* + * Set the smoothing weights, in preparation for the actual fit. + */ + h = 0.0; m = 0; + for (i=0; i<lfd->n; i++) + { wght(des,i) = weightmm(des->cf, des->wd[i],d_xi(des,i),gam); + if (wght(des,i)>0) + { if (dist(des,i)>h) h = dist(des,i); + des->ind[m] = i; + m++; + } + } + des->n = m; + return(h); +} +/* + * Copyright 1996-2006 Catherine Loader. + */ +/* + * + * Defines the weight functions and related quantities used + * in LOCFIT. + */ + +#include "locf.h" + +/* + * convert kernel and kernel type strings to numeric codes. + */ +#define NWFUNS 13 +static char *wfuns[NWFUNS] = { + "rectangular", "epanechnikov", "bisquare", "tricube", + "triweight", "gaussian", "triangular", "ququ", + "6cub", "minimax", "exponential", "maclean", "parametric" }; +static int wvals[NWFUNS] = { WRECT, WEPAN, WBISQ, WTCUB, + WTRWT, WGAUS, WTRIA, WQUQU, W6CUB, WMINM, WEXPL, WMACL, WPARM }; +int lfkernel(char *z) +{ return(pmatch(z, wfuns, wvals, NWFUNS, WTCUB)); +} + +#define NKTYPE 5 +static char *ktype[NKTYPE] = { "spherical", "product", "center", "lm", "zeon" }; +static int kvals[NKTYPE] = { KSPH, KPROD, KCE, KLM, KZEON }; +int lfketype(char *z) +{ return(pmatch(z, ktype, kvals, NKTYPE, KSPH)); +} + +/* The weight functions themselves. Used everywhere. */ +double W(u,ker) +double u; +int ker; +{ u = fabs(u); + switch(ker) + { case WRECT: return((u>1) ? 0.0 : 1.0); + case WEPAN: return((u>1) ? 0.0 : 1-u*u); + case WBISQ: if (u>1) return(0.0); + u = 1-u*u; return(u*u); + case WTCUB: if (u>1) return(0.0); + u = 1-u*u*u; return(u*u*u); + case WTRWT: if (u>1) return(0.0); + u = 1-u*u; return(u*u*u); + case WQUQU: if (u>1) return(0.0); + u = 1-u*u; return(u*u*u*u); + case WTRIA: if (u>1) return(0.0); + return(1-u); + case W6CUB: if (u>1) return(0.0); + u = 1-u*u*u; u = u*u*u; return(u*u); + case WGAUS: return(exp(-SQR(GFACT*u)/2.0)); + case WEXPL: return(exp(-EFACT*u)); + case WMACL: return(1/((u+1.0e-100)*(u+1.0e-100))); + case WMINM: LERR(("WMINM in W")); + return(0.0); + case WPARM: return(1.0); + } + LERR(("W(): Unknown kernel %d\n",ker)); + return(1.0); +} + +int iscompact(ker) +int ker; +{ if ((ker==WEXPL) | (ker==WGAUS) | (ker==WMACL) | (ker==WPARM)) return(0); + return(1); +} + +double weightprod(lfd,u,h,ker) +lfdata *lfd; +double *u, h; +int ker; +{ int i; + double sc, w; + w = 1.0; + for (i=0; i<lfd->d; i++) + { sc = lfd->sca[i]; + switch(lfd->sty[i]) + { case STLEFT: + if (u[i]>0) return(0.0); + w *= W(-u[i]/(h*sc),ker); + break; + case STRIGH: + if (u[i]<0) return(0.0); + w *= W(u[i]/(h*sc),ker); + break; + case STANGL: + w *= W(2*fabs(sin(u[i]/(2*sc)))/h,ker); + break; + case STCPAR: + break; + default: + w *= W(fabs(u[i])/(h*sc),ker); + } + if (w==0.0) return(w); + } + return(w); +} + +double weightsph(lfd,u,h,ker, hasdi,di) +lfdata *lfd; +double *u, h, di; +int ker, hasdi; +{ int i; + + if (!hasdi) di = rho(u,lfd->sca,lfd->d,KSPH,lfd->sty); + + for (i=0; i<lfd->d; i++) + { if ((lfd->sty[i]==STLEFT) && (u[i]>0.0)) return(0.0); + if ((lfd->sty[i]==STRIGH) && (u[i]<0.0)) return(0.0); + } + if (h==0) return((di==0.0) ? 1.0 : 0.0); + + return(W(di/h,ker)); +} + +double weight(lfd,sp,x,t,h, hasdi,di) +lfdata *lfd; +smpar *sp; +double *x, *t, h, di; +int hasdi; +{ double u[MXDIM]; + int i; + for (i=0; i<lfd->d; i++) u[i] = (t==NULL) ? x[i] : x[i]-t[i]; + switch(kt(sp)) + { case KPROD: return(weightprod(lfd,u,h,ker(sp))); + case KSPH: return(weightsph(lfd,u,h,ker(sp), hasdi,di)); + } + LERR(("weight: unknown kernel type %d",kt(sp))); + return(1.0); +} + +double sgn(x) +double x; +{ if (x>0) return(1.0); + if (x<0) return(-1.0); + return(0.0); +} + +double WdW(u,ker) /* W'(u)/W(u) */ +double u; +int ker; +{ double eps=1.0e-10; + if (ker==WGAUS) return(-GFACT*GFACT*u); + if (ker==WPARM) return(0.0); + if (fabs(u)>=1) return(0.0); + switch(ker) + { case WRECT: return(0.0); + case WTRIA: return(-sgn(u)/(1-fabs(u)+eps)); + case WEPAN: return(-2*u/(1-u*u+eps)); + case WBISQ: return(-4*u/(1-u*u+eps)); + case WTRWT: return(-6*u/(1-u*u+eps)); + case WTCUB: return(-9*sgn(u)*u*u/(1-u*u*fabs(u)+eps)); + case WEXPL: return((u>0) ? -EFACT : EFACT); + } + LERR(("WdW: invalid kernel")); + return(0.0); +} + +/* deriv. weights .. spherical, product etc + u, sc, sty needed only in relevant direction + Acutally, returns (d/dx W(||x||/h) ) / W(.) +*/ +double weightd(u,sc,d,ker,kt,h,sty,di) +double u, sc, h, di; +int d, ker, kt, sty; +{ if (sty==STANGL) + { if (kt==KPROD) + return(-WdW(2*sin(u/(2*sc)),ker)*cos(u/(2*sc))/(h*sc)); + if (di==0.0) return(0.0); + return(-WdW(di/h,ker)*sin(u/sc)/(h*sc*di)); + } + if (sty==STCPAR) return(0.0); + if (kt==KPROD) + return(-WdW(u/(h*sc),ker)/(h*sc)); + if (di==0.0) return(0.0); + return(-WdW(di/h,ker)*u/(h*di*sc*sc)); +} + +double weightdd(u,sc,d,ker,kt,h,sty,di,i0,i1) +double *u, *sc, h, di; +int d, ker, kt, i0, i1, *sty; +{ double w; + w = 1; + if (kt==KPROD) + { + w = WdW(u[i0]/(h*sc[i0]),ker)*WdW(u[i1]/(h*sc[i1]),ker)/(h*h*sc[i0]*sc[i1]); + } + return(0.0); +} + +/* Derivatives W'(u)/u. + Used in simult. conf. band computations, + and kernel density bandwidth selectors. */ +double Wd(u,ker) +double u; +int ker; +{ double v; + if (ker==WGAUS) return(-SQR(GFACT)*exp(-SQR(GFACT*u)/2)); + if (ker==WPARM) return(0.0); + if (fabs(u)>1) return(0.0); + switch(ker) + { case WEPAN: return(-2.0); + case WBISQ: return(-4*(1-u*u)); + case WTCUB: v = 1-u*u*u; + return(-9*v*v*u); + case WTRWT: v = 1-u*u; + return(-6*v*v); + default: LERR(("Invalid kernel %d in Wd",ker)); + } + return(0.0); +} + +/* Second derivatives W''(u)-W'(u)/u. + used in simult. conf. band computations in >1 dimension. */ +double Wdd(u,ker) +double u; +int ker; +{ double v; + if (ker==WGAUS) return(SQR(u*GFACT*GFACT)*exp(-SQR(u*GFACT)/2)); + if (ker==WPARM) return(0.0); + if (u>1) return(0.0); + switch(ker) + { case WBISQ: return(12*u*u); + case WTCUB: v = 1-u*u*u; + return(-9*u*v*v+54*u*u*u*u*v); + case WTRWT: return(24*u*u*(1-u*u)); + default: LERR(("Invalid kernel %d in Wdd",ker)); + } + return(0.0); +} + +/* int u1^j1..ud^jd W(u) du. + Used for local log-linear density estimation. + Assume all j_i are even. + Also in some bandwidth selection. +*/ +double wint(d,j,nj,ker) +int d, *j, nj, ker; +{ double I, z; + int k, dj; + dj = d; + for (k=0; k<nj; k++) dj += j[k]; + switch(ker) /* int_0^1 u^(dj-1) W(u)du */ + { case WRECT: I = 1.0/dj; break; + case WEPAN: I = 2.0/(dj*(dj+2)); break; + case WBISQ: I = 8.0/(dj*(dj+2)*(dj+4)); break; + case WTCUB: I = 162.0/(dj*(dj+3)*(dj+6)*(dj+9)); break; + case WTRWT: I = 48.0/(dj*(dj+2)*(dj+4)*(dj+6)); break; + case WTRIA: I = 1.0/(dj*(dj+1)); break; + case WQUQU: I = 384.0/(dj*(dj+2)*(dj+4)*(dj+6)*(dj+8)); break; + case W6CUB: I = 524880.0/(dj*(dj+3)*(dj+6)*(dj+9)*(dj+12)*(dj+15)*(dj+18)); break; + case WGAUS: switch(d) + { case 1: I = S2PI/GFACT; break; + case 2: I = 2*PI/(GFACT*GFACT); break; + default: I = exp(d*log(S2PI/GFACT)); /* for nj=0 */ + } + for (k=0; k<nj; k++) /* deliberate drop */ + switch(j[k]) + { case 4: I *= 3.0/(GFACT*GFACT); + case 2: I /= GFACT*GFACT; + } + return(I); + case WEXPL: I = factorial(dj-1)/ipower(EFACT,dj); break; + default: LERR(("Unknown kernel %d in exacint",ker)); + } + if ((d==1) && (nj==0)) return(2*I); /* common case quick */ + z = (d-nj)*LOGPI/2-mut_lgammai(dj); + for (k=0; k<nj; k++) z += mut_lgammai(j[k]+1); + return(2*I*exp(z)); +} + +/* taylor series expansion of weight function around x. + 0 and 1 are common arguments, so are worth programming + as special cases. + Used in density estimation. +*/ +int wtaylor(f,x,ker) +double *f, x; +int ker; +{ double v; + switch(ker) + { case WRECT: + f[0] = 1.0; + return(1); + case WEPAN: + f[0] = 1-x*x; f[1] = -2*x; f[2] = -1; + return(3); + case WBISQ: + v = 1-x*x; + f[0] = v*v; f[1] = -4*x*v; f[2] = 4-6*v; + f[3] = 4*x; f[4] = 1; + return(5); + case WTCUB: + if (x==1.0) + { f[0] = f[1] = f[2] = 0; f[3] = -27; f[4] = -81; f[5] = -108; + f[6] = -81; f[7] = -36; f[8] = -9; f[9] = -1; return(10); } + if (x==0.0) + { f[1] = f[2] = f[4] = f[5] = f[7] = f[8] = 0; + f[0] = 1; f[3] = -3; f[6] = 3; f[9] = -1; return(10); } + v = 1-x*x*x; + f[0] = v*v*v; f[1] = -9*v*v*x*x; f[2] = x*v*(27-36*v); + f[3] = -27+v*(108-84*v); f[4] = -3*x*x*(27-42*v); + f[5] = x*(-108+126*v); f[6] = -81+84*v; + f[7] = -36*x*x; f[8] = -9*x; f[9] = -1; + return(10); + case WTRWT: + v = 1-x*x; + f[0] = v*v*v; f[1] = -6*x*v*v; f[2] = v*(12-15*v); + f[3] = x*(20*v-8); f[4] = 15*v-12; f[5] = -6; f[6] = -1; + return(7); + case WTRIA: + f[0] = 1-x; f[1] = -1; + return(2); + case WQUQU: + v = 1-x*x; + f[0] = v*v*v*v; f[1] = -8*x*v*v*v; f[2] = v*v*(24-28*v); + f[3] = v*x*(56*v-32); f[4] = (70*v-80)*v+16; f[5] = x*(32-56*v); + f[6] = 24-28*v; f[7] = 8*x; f[8] = 1; + return(9); + case W6CUB: + v = 1-x*x*x; + f[0] = v*v*v*v*v*v; + f[1] = -18*x*x*v*v*v*v*v; + f[2] = x*v*v*v*v*(135-153*v); + f[3] = v*v*v*(-540+v*(1350-816*v)); + f[4] = x*x*v*v*(1215-v*(4050-v*3060)); + f[5] = x*v*(-1458+v*(9234+v*(-16254+v*8568))); + f[6] = 729-v*(10206-v*(35154-v*(44226-v*18564))); + f[7] = x*x*(4374-v*(30132-v*(56862-v*31824))); + f[8] = x*(12393-v*(61479-v*(92664-v*43758))); + f[9] = 21870-v*(89100-v*(115830-v*48620)); + f[10]= x*x*(26730-v*(69498-v*43758)); + f[11]= x*(23814-v*(55458-v*31824)); + f[12]= 15849-v*(34398-v*18564); + f[13]= x*x*(7938-8568*v); + f[14]= x*(2970-3060*v); + f[15]= 810-816*v; + f[16]= 153*x*x; + f[17]= 18*x; + f[18]= 1; + return(19); + } + LERR(("Invalid kernel %d in wtaylor",ker)); + return(0); +} + +/* convolution int W(x)W(x+v)dx. + used in kde bandwidth selection. +*/ +double Wconv(v,ker) +double v; +int ker; +{ double v2; + switch(ker) + { case WGAUS: return(SQRPI/GFACT*exp(-SQR(GFACT*v)/4)); + case WRECT: + v = fabs(v); + if (v>2) return(0.0); + return(2-v); + case WEPAN: + v = fabs(v); + if (v>2) return(0.0); + return((2-v)*(16+v*(8-v*(16-v*(2+v))))/30); + case WBISQ: + v = fabs(v); + if (v>2) return(0.0); + v2 = 2-v; + return(v2*v2*v2*v2*v2*(16+v*(40+v*(36+v*(10+v))))/630); + } + LERR(("Wconv not implemented for kernel %d",ker)); + return(0.0); +} + +/* derivative of Wconv. + 1/v d/dv int W(x)W(x+v)dx + used in kde bandwidth selection. +*/ +double Wconv1(v,ker) +double v; +int ker; +{ double v2; + v = fabs(v); + switch(ker) + { case WGAUS: return(-0.5*SQRPI*GFACT*exp(-SQR(GFACT*v)/4)); + case WRECT: + if (v>2) return(0.0); + return(1.0); + case WEPAN: + if (v>2) return(0.0); + return((-16+v*(12-v*v))/6); + case WBISQ: + if (v>2) return(0.0); + v2 = 2-v; + return(-v2*v2*v2*v2*(32+v*(64+v*(24+v*3)))/210); + } + LERR(("Wconv1 not implemented for kernel %d",ker)); + return(0.0); +} + +/* 4th derivative of Wconv. + used in kde bandwidth selection (BCV, SJPI, GKK) +*/ +double Wconv4(v,ker) +double v; +int ker; +{ double gv; + switch(ker) + { case WGAUS: + gv = GFACT*v; + return(exp(-SQR(gv)/4)*GFACT*GFACT*GFACT*(12-gv*gv*(12-gv*gv))*SQRPI/16); + } + LERR(("Wconv4 not implemented for kernel %d",ker)); + return(0.0); +} + +/* 5th derivative of Wconv. + used in kde bandwidth selection (BCV method only) +*/ +double Wconv5(v,ker) /* (d/dv)^5 int W(x)W(x+v)dx */ +double v; +int ker; +{ double gv; + switch(ker) + { case WGAUS: + gv = GFACT*v; + return(-exp(-SQR(gv)/4)*GFACT*GFACT*GFACT*GFACT*gv*(60-gv*gv*(20-gv*gv))*SQRPI/32); + } + LERR(("Wconv5 not implemented for kernel %d",ker)); + return(0.0); +} + +/* 6th derivative of Wconv. + used in kde bandwidth selection (SJPI) +*/ +double Wconv6(v,ker) +double v; +int ker; +{ double gv, z; + switch(ker) + { case WGAUS: + gv = GFACT*v; + gv = gv*gv; + z = exp(-gv/4)*(-120+gv*(180-gv*(30-gv)))*0.02769459142; + gv = GFACT*GFACT; + return(z*gv*gv*GFACT); + } + LERR(("Wconv6 not implemented for kernel %d",ker)); + return(0.0); +} + +/* int W(v)^2 dv / (int v^2 W(v) dv)^2 + used in some bandwidth selectors +*/ +double Wikk(ker,deg) +int ker, deg; +{ switch(deg) + { case 0: + case 1: /* int W(v)^2 dv / (int v^2 W(v) dv)^2 */ + switch(ker) + { case WRECT: return(4.5); + case WEPAN: return(15.0); + case WBISQ: return(35.0); + case WGAUS: return(0.2820947918*GFACT*GFACT*GFACT*GFACT*GFACT); + case WTCUB: return(34.152111046847892); /* 59049 / 1729 */ + case WTRWT: return(66.083916083916080); /* 9450/143 */ + } + case 2: + case 3: /* 4!^2/8*int(W1^2)/int(v^4W1)^2 + W1=W*(n4-v^2n2)/(n0n4-n2n2) */ + switch(ker) + { case WRECT: return(11025.0); + case WEPAN: return(39690.0); + case WBISQ: return(110346.9231); + case WGAUS: return(14527.43412); + case WTCUB: return(126500.5904); + case WTRWT: return(254371.7647); + } + } + LERR(("Wikk not implemented for kernel %d, deg %d",ker,deg)); + return(0.0); +}