new; cls;
label={"A","B","C","D","E","F"};
a={14,23,2.2,17,1.6,1.10};
b={15,25,2.5,23,3.2,1.15};
c={19,33,3.1,35,3.3,1.17};
times=1000;       


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


/*
**  sen2md.txt - Calculation of Sensitivities from each variable by contrapositive simulation.
** (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 sensitivities from each variable assuming triang(a,b,c)
**              in terms of mean deviation by the contrapositive simulation of sen1md.
**
**  Format:     call sen2md(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:     dsadj    n x 1 vector, normalized change of mean deviations
**
**
**  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.
**
**              This contrapositive method is also valid and accurate when we ignore correlation among variables
**              or assume zezo correlation.
**
*/

proc sen2md(label,a,b,c,times);
   local n,m,i,j,sd0,sd,results,mm,ds,dsadj,mask;
   n=rows(a);
   m=zeros(n,times);
   j=1;
   do while j<=times;
      i=1;
      do while i<=n;
         m[i,j]=(a[i]+b[i]+c[i])/3;
         i=i+1;
      endo;
      j=j+1;
   endo;
   sd0=0;
   sd=zeros(n,1);
   i=1;
   do while i<=n;
      mm=m;
      mm[i,.]=triang2(a[i],b[i],c[i],times)';
      results=model(mm);
      sd[i]=md(results');
      i=i+1;
   endo;
   ds=-(sd0-sd);                                       /* now, negative of them */
   dsadj=ds/sumc(ds);
/* display the results */
   mask={0 1 1 1};
   print "                     m.d.          d(m.d.)    normalized d(m.d.) %";
   call printfmt((miss(0,0)|label)~(sd0|sd)~(miss(0,0)|ds)~(miss(0,0)|dsadj*100),mask);
   print/lz "                            sum:" sumc(ds);
   print;
/* graph */
   library pgraph;
   graphset;
   pqgwin auto;
   ylabel("% Sensitivity");
   _pdate="";
   _protate=1;
   bar(0,dsadj*100);
   retp(dsadj);
endp;

proc md(x);
   local n,xbar;
   n=rows(x);
   xbar=meanc(x);
   retp( sumc(abs(x-xbar'))/n );
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;