/* Bezier Curve */
/* 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) */
call beziercurve(y,x,allpoints);




proc beziercurve(y,x,allpoints);
   local z,n,xx,yy,t,i;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   n=rows(x);
/* t step (0<=1<=1) */
   t=seqa(0,1/(allpoints-1),allpoints); t[allpoints]=1;
   {yy,xx}=bezier(y,x,t);
/* graph */
   library pgraph;
   graphset;
   pqgwin auto;
   begwind;
   window(1,1,0);
   scale(-0.05*(maxc(x)-minc(x))+minc(x)|maxc(x)+0.05*(maxc(x)-minc(x)),-0.25*(maxc(y)-minc(y))+minc(y)|maxc(y)+0.25*(maxc(y)-minc(y)));
   setwind(1);
      title("Bezier Curve");
      xy(xx,yy);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(x,y);
   endwind;
   retp(xx~yy);
endp;

proc(2)=bezier(y,x,t);
   local z,n,yi,xi,i;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   n=rows(x);
/* main */
   n=n-1;
   yi=0; xi=0;
   i=0;
   do while i<=n;
      yi=yi+nCr(n,i)*(1-t)^(n-i).*t^i*y[i+1];
      xi=xi+nCr(n,i)*(1-t)^(n-i).*t^i*x[i+1];
      i=i+1;
   endo;
   retp(yi,xi);
endp;

proc nCr(n,r);
   if r==0 or n==r;
      retp(1);
   endif;
   retp( exp(sumc(ln(seqa(n,-1,n)))-sumc(ln(seqa(r,-1,r)))-sumc(ln(seqa(n-r,-1,n-r)))) );
endp;