proc iml ;
     
* Simulation Parameters ;
    
    n=100000;  * Number of simulated data sets ;
    c=20;       * Number of components ;
    h=150;      * Number of Hours;
    t=15;       * Number need to survive ;
    lambda=300; * Average lifetime;
    
    count=0;    * Number of times at least t components survive ;
*   results=j(n,1,.); *Results for each data set;
    do run=1 to n;  * Data set loop ;
        time=j(c,1,.);
        do comp=1 to c;
            time[comp]=lambda*ranexp(-1); * Generate lifetime for each component ;
            end;
        Surv=(time>=h);                * Component survives h hours ;
        Nsurv=sum(Surv);               * Number of components that survive ;
        count=count+(Nsurv>=t);        * At least t components have survived;
*       results[run]=(Nsurv>=t);
        end;
    Prob=count/n;  * Probability of at least t components surviving;
    Print Prob[label="Estimated Probability" format=5.4];
*   Exact Probability ;
    p1=exp(-h/lambda);
    p2=0; * Initialize ;
    do i=t to c;
        p2=p2+comb(c,i)*p1**i*(1-p1)**(c-i);
        end;
    Print P2[label="Exact Probability" format=5.4];
    quit;
run;