Estimates parameters for beta distribution from mean and sd
Source:    Publish Time: 2012-08-24 08:10   1326 Views   Size:  16px  14px  12px
Author: Xuanqian Xie  *1. It was assumed that probabilities and utility are based on beta distribution. We know

Author: Xuanqian Xie 

*1. It was assumed that probabilities and utility are based on beta distribution.

We know mean and sd. For the PSA purpose, we need know the parameters a and b of the beta distribution, which are n nonlinear equations in n unknowns.

Beta distribution (a, b). Mean= a/(a+b)  Sd=((a*b)/(a+b)2(a+b+1))1/2.

I use the SOLVE statement in PROC MODEL to obtain the unknowns, which applies Newton's method for solving equations.

The theory of Newton's method can be found at

http://library.thinkquest.org/C006002/Pages/Newtons_Method_for_Solving_Equations.htm .

We also need provide initial values of equations.

The initial values must be appropriate.

Otherwise, the equations can not be solved for beyond extreme value, and so on.

Newton's method of SAS program can provide one solution of some possible solutions.*;

 

%macro Newton (range=&ra, No=10000, out=Prob1);

*Import data;

%let pa=%str(C:\Documents and Settings\sxie\My Documents\Maze\MAZE-model\Inputs for Markov model-new.xls);

proc import

    datafile="&pa"

    out=aa

    DBMS=EXCEL2000 REPLACE;

    getnames=yes;

    range="&range";

run;

 

* Manage data;

Data bb;

    set aa;

      b=SD_PSA_*SD_PSA_;

    x1=prob*&no;

    x2=(1-prob)*&no;

      keep name cycle_No prob b x1 x2;

      where SD_PSA_ ^=.;

run;

 

*Calculate parameters of Beta distribution;

proc model data=bb;

    eq.one= (x1/(x1+x2))-prob;

    eq.two= (x1*x2)/((x1+x2)*(x1+x2)*(x1+x2+1))-b;

    solve x1 x2/itprint out=cc outpredict;

quit;

run;

 

*Combine data set for output;

proc sql;

   create table dd as

    select distinct name, a.cycle_No, a.Prob, x1 as a, x2 as b

        from aa as a

          left join

         cc as c

        on a.Prob = c.Prob

   where name ^=""

   order by a.name, a.cycle_No ;

quit;

 

*output data set to excel;

proc export data=dd

   outfile="&pa"

   DBMS=Excel2000 replace;

   sheet="&out";

run;

 

*Delete processing files;

proc datasets lib=work;

    delete aa bb cc dd ;

run;

quit;

 

%mend;

 

%let ra=%str(Prob$a3:g72);

%Newton (range=&ra, No=10000, out=Prob1);

 

%let ra721=%str(Prob$a3:g16);

%Newton (range=&ra721, No=10000, out=Prob721);

 

%let ut=%str(utility$a89:d106);

%Newton (range=&ut, No=100, out=utility1);