topic Re: Rand() in SAS Programming

Rand()

Azeddine — Wed, 14 Sep 2016 19:42:42 GMT

Hello,

I use Proc severity and i have parameter for GPD distribution,

I would like to use them to simulate an sampling of GPD like a "rand("Weib", 0.25,2.1)"

My parameters are for theta =126998.5 and for xi =0.31954.

Re: Rand()

Rick_SAS — Thu, 15 Sep 2016 14:08:49 GMT

Reference: p 113 of Wicklin (2013) Simulating Data with SAS

The discussion begins with "you can use the inverse CDF algorithm to simulate random variates." It ends with this DATA step code:

U = rand("Uniform");

X = theta - sigma/alpha * (U**alpha-1);

The three-parameter distribution for the above uses the parameterization documented in PROC UNIVARIATE: theta=threshold, sigma=scale, and alpha=shape.

Re: Rand()

Azeddine — Tue, 20 Sep 2016 07:51:48 GMT

Thank you for your answer, in my case it doesn't work,

So i try with a analytic inversion of function:

The initial Function of GPD is defined here in proc severity:

http://support.sas.com/documentation/cdl/en/etsug/63348/HTML/default/viewer.htm#etsug_severity_sect015.htm

with x=(x'-Treshold)

I use this formula :

X'=&Treshold+(((1-rand('uniform'))**(-&&Gpd_xi))-1)*(&&Gpd_theta/&&Gpd_xi)

Re: Rand()

Rick_SAS — Tue, 20 Sep 2016 09:54:01 GMT

You say it doesn't work. Please show us the SAS code you are using.

Re: Rand()

Azeddine — Tue, 20 Sep 2016 10:08:55 GMT

To be honnest it gave me result but the result are different between resimulation with proc univariate and your resimulation.

Re: Rand()

Rick_SAS — Tue, 20 Sep 2016 10:18:34 GMT

If you shows that DATA step and PROC UNIVARIATE call, we can help diagnose the problem.

Re: Rand()

Rick_SAS — Tue, 20 Sep 2016 13:18:19 GMT

Here is simulation code and PROC UINIVARIATE code that shows agreement between the parameters and the parameter estimates.

/* simulation of generalized pareto distribution */
data GenPareto(keep= X);
alpha = 1/4;     /* shape */   
sigma = 10;      /* scale */
theta = 1.5;     /* threshold */
call streaminit(1);
/* support of x is theta + (0, sigma/alpha] = theta + (0,40] */
do i = 1 to 400;
   U = rand("Uniform");
   X = theta - sigma/alpha * (U**alpha - 1);
   output;
end;
run;

proc univariate data=GenPareto;
histogram x / pareto(threshold=1.5); /* parameters are alpha=0.25 sigma=10 */
ods select histogram;
run;