Rand()

Accepted Solution Solved
Reply
Occasional Contributor
Posts: 18
Accepted Solution

Rand()

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.

 

 


Accepted Solutions
Solution
‎09-20-2016 06:19 AM
SAS Super FREQ
Posts: 3,637

Re: Rand()

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

View solution in original post


All Replies
SAS Super FREQ
Posts: 3,637

Re: Rand()

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.

Occasional Contributor
Posts: 18

Re: Rand()

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_sect0...

with x=(x'-Treshold)

 

I use this formula :

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

 

SAS Super FREQ
Posts: 3,637

Re: Rand()

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

 

Occasional Contributor
Posts: 18

Re: Rand()

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

 

 

Solution
‎09-20-2016 06:19 AM
SAS Super FREQ
Posts: 3,637

Re: Rand()

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

SAS Super FREQ
Posts: 3,637

Re: Rand()

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;
☑ This topic is solved.

Need further help from the community? Please ask a new question.

Discussion stats
  • 6 replies
  • 723 views
  • 0 likes
  • 2 in conversation