new; cls;
label={"A","B","C","D","E","F","G"};
a={12,13,4,3,4,11,5};
b={17,23,7.5,7.5,5.5,14,8.5};
c={22,33,11,12,7,17,12};
times=1000;

fn model(x)=x[1,.]+x[2,.]+x[3,.]+x[4,.]+x[5,.]+x[6,.]+x[7,.];        /*  model setting */
call costmodel1(label,a,b,c,times);



/*
**  contrib1.txt - Calculation of Contributions from each variable by Rank Correlation.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Reference:  David Vose
**              "Risk Analysis: A Quantitative Guide" 2nd Edition, Ch.15, John Wiley & Sons 2000.
**
**
**  Purpose:    Computes contributions from each variable by rank correlation coefficients
**              assuming triang(a,b,c) for each variable.
**
**  Format:     call costmodel1(label,a,b,c,times);
**
**  Input:      label    n x 1 vector, labels;
**
**              a        n x 1 vector, minimums of 1st,2nd,3rd,..., nth variable
**
**              b        n x 1 vector, most likely values of 1st,...nth variable
**
**              c        n x 1 vector, maximums of 1st,2nd,3rd,..., nth variable
**
**              times    scalar, number of simulations
**
**  Output:     normrho  n x 1 vector, normalized rank correlation coefficients(contributions)
**
**
**  Notice:     You need to set the model you are trying to evaluate before calling this procedure
**              by vector notations, i.e. you must use '.*' and './' for multiplication and division.
**
**              Normalized range is calculated with signs which would be the precise way when there
**              are some negative rank correlation coefficients. If you don't like that, replace
**              range1=range.*((-1)^(rho.<0)); to range1=range; .
*/

proc costmodel1(label,a,b,c,times);
   local n,dist,i,sum,rho,range,range1,normrange,normrho,mask,fmt;
   n=rows(a);
   dist=zeros(n,times);
   i=1;
   do while i<=n;
      dist[i,.]=triang2(a[i],b[i],c[i],times)';
      i=i+1;
   endo;
   sum=model(dist);
   sum=sum';
screen off;
   rho=zeros(n,1);
   i=1;
   do while i<=n;
      rho[i]=speaman(sum,dist[i,.]');             /* Or, you could use speaman2. */
      i=i+1;
   endo;
screen on;
   range=c-a;
   range1=range.*((-1)^(rho.<0));
   normrange=range1/sumc(range1);
   normrho=rho/sumc(rho);
/* display the results */
   print "                                                                                  Normalized";
   print "   Component     Minimum Most likely     Maximum       Range   Rank corr.      Range   Rank corr.";
   mask={0 1 1 1 1 1 1 1};
   fmt="*.*lg"~12~6;
   call printfm(label~a~b~c~range~rho~normrange*100~normrho*100,mask,fmt);
   format/rz 11,6;
   print "                                             SUM:" sumc(range) sumc(rho) sumc(normrange*100) sumc(normrho*100);
   print;
   format /mb1 /ros 16,8;
/* graph */
   library pgraph;
   graphset;
   pqgwin auto;
   ylabel("% Normalized Rank Correlation Coefficients");
   _pdate="";
   _protate=1;
   bar(0,normrho*100);
   retp(normrho);
endp;

proc speaman(y,x);
   local n,u,v,rho;
   n=rows(x);
   u=ranksame(x);
   v=ranksame(y);
   rho=1-6*sumc((u-v)^2)/(n*(n^2-1));      /* (3.25) */
   print/rz x~y~u~v~(u-v)^2;
   print;
   print/rz "  n=" n;
   print/rz "rho=" rho;
   print;
   retp(rho);
endp;

proc speaman2(y,x);
   local n,u,v,SSuv,SSuu,SSvv,rho;
   n=rows(x);
   u=ranksame(x);
   v=ranksame(y);
   SSuv=sumc((u-meanc(u)).*(v-meanc(v)));
   SSuu=sumc((u-meanc(u))^2);
   SSvv=sumc((v-meanc(v))^2);
   rho=SSuv/sqrt(SSuu*SSvv);               /* (3.26) */
   print/rz x~y~u~v~(u-v)^2;
   print;
   print/rz "  n=" n;
   print/rz "rho=" rho;
   print;
   retp(rho);
endp;

proc ranksame(x);
   local n,y,i,j,k,index;
   n=rows(x);
   y=zeros(n,1);
   i=1;
   do while i<=n;
      index=sumc(x.==minc(x));
      k=0;
      j=0;
      do while j<=index-1;
         k=k+(i+j);
         j=j+1;
      endo;
      k=k/index;
      y=y+k*(x.==minc(x));
      x=x+(x.==minc(x))*1e256;
      i=i+index;
   endo;
   retp(y);
endp;





proc triang(a,b,c,n);
   local x,i;
   if b==c;
      x=rndu(n,1);
      x=a+sqrt(x)*(b-a);
   elseif a==b;
      x=rndu(n,1);
      x=b+(1-sqrt(x))*(c-b);
   else;
      x=rndu(n,1);
      i=1;
      do while i<=n;
         if rndu(1,1)<(b-a)/(c-a);
            x[i]=a+sqrt(x[i])*(b-a);
         else;
            x[i]=b+(1-sqrt(x[i]))*(c-b);
         endif;
         i=i+1;
      endo;
   endif;
   retp(x);
endp;

proc triang2(a,b,c,n);
   local n1,n2,d,U,i,x,x1,x2,xx,index;
   if b==c;
      d=1/n;
      U=zeros(n,1);
      i=1;
      do while i<=n;
         U[i]=(d*i-d*(i-1))*rndu(1,1)+d*(i-1);
         i=i+1;
      endo;
      x=a+sqrt(U)*(b-a);
   elseif a==b;
      d=1/n;
      U=zeros(n,1);
      i=1;
      do while i<=n;
         U[i]=(d*i-d*(i-1))*rndu(1,1)+d*(i-1);
         i=i+1;
      endo;
      x=b+(1-sqrt(U))*(c-b);
   else;  
      n1=round(n*(b-a)/(c-a));
      n2=n-n1;
   /* a to b */
      d=1/n1;
      U=zeros(n1,1);
      i=1;
      do while i<=n1;
         U[i]=(d*i-d*(i-1))*rndu(1,1)+d*(i-1);
         i=i+1;
      endo;
      x1=a+sqrt(U)*(b-a);
   /* b to c */
      d=1/n2;
      U=zeros(n2,1);
      i=1;
      do while i<=n2;
         U[i]=(d*i-d*(i-1))*rndu(1,1)+d*(i-1);
         i=i+1;
      endo;
      x2=b+(1-sqrt(U))*(c-b);
      x=x1|x2;
   endif;
/* randomize the order  */
   xx=zeros(n,1);
   i=1;
   do while i<=n-1;
      index=ceil((n-(i-1))*rndu(1,1));
      xx[i]=x[index];
      if index==1;
         x=x[2:rows(x)];
      elseif index==rows(x);
         x=x[1:rows(x)-1];
      else;
         x=x[1:(index-1) (index+1):rows(x)];
      endif;
      i=i+1;
   endo;
   xx[n]=x;
   x=xx;
   retp(x);
endp;