/* NURBS eliminating some points */
/* Data from Holton(2003)"Value-at-Risk" */
new; cls;
let data[10,2]=
1.1  2.14
1.4  2.60
2.5  1.15
2.7  1.19
3.2  1.88
3.6  1.55
4.1  2.65
4.3  3.80
4.5  4.46
4.9  6.35
;
x=data[.,1]; y=data[.,2];
allpoints=100;                               /* # of all points between min(x) and max(x) */
k=3;                                         /* # of order */
w=ones(rows(x),1);                           /* n x 1 weight vector */
eliminate={1,2};
call nurbs1(y,x,w,k,allpoints,eliminate);


proc nurbs1(y,x,w,k,allpoints,eliminate);
   local z,n,x0,i,j,xx,yy,d,x1,xx1,b1,y1,x_1,xx_1,b_1,y_1,min,max,sumw,eindex,xall,yall;
/* sort them out in terms of x */
   z=x~y~w;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   w=z[.,3];
/* elimination */
   if maxc(eliminate)>rows(x);
      errorlog "ERROR: Out of range.";
      retp(".");
   elseif eliminate==0 or eliminate==zeros(rows(eliminate),1);
      errorlog "No elimination:";
      xall=x;
      yall=y;
      w=ones(k+1,1)|w|ones(k+1,1);
   else;
      xall=x;
      yall=y;
      eindex=zeros(rows(x),1);
      eindex[eliminate]=ones(rows(eliminate),1);
      x=delif(x,eindex);
      y=delif(y,eindex);
      w=delif(w,eindex);
      w=ones(k+1,1)|w|ones(k+1,1);
   endif;
/* extend points */
   n=rows(x);
   d=x[n-k:n]-x[n-k-1];
   x1=x[n]+d;
   xx1=x1^seqa(0,1,k+1)';
   b1=inv(xx1)*y[n-k:n];
   y1=xx1*b1;
   x=x|x1; y=y|y1;
   d=x[2:k+2]-x[1];
   x_1=x[1]-d;
   xx_1=x_1^seqa(0,1,k+1)';
   b_1=inv(xx_1)*y[1:k+1];
   y_1=xx_1*b_1;
   x=rev(x_1)|x; y=rev(y_1)|y;
/* main */
   min=(x[k+2]+x[k+3])/2; max=(x[k+1+n]+x[k+2+n])/2;
   x0=min+(max-min)/(allpoints-1)*seqa(0,1,allpoints);
   xx=miss(zeros(allpoints,1),0); yy=miss(zeros(allpoints,1),0);
   j=1;
   do while j<=allpoints;
      xx[j]=0; yy[j]=0; sumw=0;
      i=1;
      do while i<=rows(x)-k+1;
         xx[j]=xx[j]+w[i]*x[i]*dB(i+1-(k+1)/2,k,x,x0[j]);
         yy[j]=yy[j]+w[i]*y[i]*dB(i+1-(k+1)/2,k,x,x0[j]);
         sumw=sumw+w[i]*dB(i+1-(k+1)/2,k,x,x0[j]);
         i=i+1;
      endo;
      xx[j]=xx[j]/sumw;
      yy[j]=yy[j]/sumw;
      j=j+1;
   endo;
/* graph */
   library pgraph;
   graphset;
   pqgwin auto;
   begwind;
   window(1,1,0);
   scale(-0.05*(maxc(xall)-minc(xall))+minc(xall)|maxc(xall)+0.05*(maxc(xall)-minc(xall)),-0.25*(maxc(yall)-minc(yall))+minc(yall)|maxc(yall)+0.25*(maxc(yall)-minc(yall)));
   setwind(1);
      title("NURBS");
      xy(xx,yy);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(xall,yall);
   endwind;
   retp(xx~yy);
endp;

proc dB(i,k,t,t0);
   local B,ii,kk,index;
   if rows(t)-1<k or i<=0 or i>=rows(t);
      retp(0);
   endif;
   B=miss(zeros(k,k),0);
   ii=1;
   do while ii<=k;
      index=i+(ii-1);
      if t[index]<=t0 and t0<t[index+1];
         B[ii,1]=1;
      else;
         B[ii,1]=0;
      endif;
      ii=ii+1;
   endo;
   kk=2;
   do while kk<=k;
      ii=1;
      do while ii<=k-(kk-1);
         index=i+(ii-1);
         B[ii,kk]=(t0-t[index])/(t[index+kk-1]-t[index])*B[ii,kk-1]+(t[index+kk]-t0)/(t[index+kk]-t[index+1])*B[ii+1,kk-1];
         ii=ii+1;
      endo;
      kk=kk+1;
   endo;
   retp(B[1,k]);
endp;