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 costmodel1p(label,a,b,c,times);



/*
**  contrib1p.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 PERT(a,b,c) instead for each variable.
**
**  Format:     call costmodel1p(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 costmodel1p(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,.]=PERT(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 PERT(a,b,c,n);
   local mu,alpha1,alpha2,x;
   mu=(a+4*b+c)/6;
   alpha1=(mu-a)*6/(c-a);               /* (2b-a-c)'s are cancelled out to be 6. */
   alpha2=alpha1*(c-mu)/(mu-a);
   x=rndbeta(n,1,alpha1,alpha2)*(c-a)+a;
   retp(x);
endp;