1.利用RANGAM,这种方法的原理是利用概率的逆运算。
我们知道卡方分布是特殊的伽马分布,伽马分布的形状参数alpha=n/2,尺度参数l=0.5时,它就是自由度为n的卡方分布,故可以通过RANGAR生成卡方随机数
例如:我们生成自由度为3的卡方分布,100个随机数
data a;
do i=1 to 100;
Ch = RanGam(121212,3/2)/(0.5);
Output;
end;
run;
2.通过CINV生成卡方随机数,这是通过累积分布函数的逆运算的原理来生成卡方随机数,只要知道累积分布函数就可以
该方法与RANGAM有重合,对于RANxxx等函数没有涉及的分布,可以考虑用这种方法,像BETA、F、t分布。
例如,以上例为例
data b;
do i=1 to 100;
ch= cinv(ranuni(121212),3);
Output;
end;
run;
3.如果累积分布函数CDF的逆计算不能进行,但累积分布函数CDF和概率密度函数PDF可以计算,可以考虑试位法。
*** Generating Chi-Squared Values ***;
data test;
/* Specify the degrees of freedom */
df=3 ;
do i=1 to 100;
/* Generate a Uniform random value */
u=ranuni(121212);
/* Null value to compare to x. */
x0=0;
/* Initialize the value of x. */
x1=1;
/* Stopping rule */
do while(abs(x1-x0)>.000005);
/* Count the number of iterations. */
count+1;
/* Save the current value of x. */
x0=x1;
/* Newton's Formula */
x1=x1 -(probchi(x1,df)-u)/pdf('CHISQUARED',x1,df);
end;
/* Use SAS's inverse Chi-Sq function */
/* for comparison. */
real=cinv(u,df);
file print;
put "estimated value= " x1 @32 count=
@45 """real value""= " real;
end;
run;
4.检验
proc univariate data=a;
var ch ;
histogram /gamma(alpha=3);
run;
Parameters for Gamma Distribution | ||
Parameter | Symbol | Estimate |
Threshold | Theta | 0 |
Scale | Sigma | 1.016927 |
Shape | Alpha | 3 |
Mean |
| 3.05078 |
Std Dev |
| 1.761369 |
Test | Statistic | p Value | ||
Kolmogorov-Smirnov | D | 0.12825562 | Pr > D | 0.003 |
Cramer-von Mises | W-Sq | 0.39848399 | Pr > W-Sq | 0.001 |
Anderson-Darling | A-Sq | 5.39904802 | Pr > A-Sq | <0.001 |
Percent | Quantile | |
Observed | Estimated | |
1.0 | 0.05840 | 0.44343 |
5.0 | 0.26045 | 0.83153 |
10.0 | 0.57196 | 1.12072 |
25.0 | 1.34261 | 1.75654 |
50.0 | 2.61212 | 2.71932 |
75.0 | 4.30397 | 3.98676 |
90.0 | 6.53834 | 5.41241 |
95.0 | 7.50933 | 6.40236 |
99.0 | 11.11379 | 8.54823 |
参考:《Pseudo-Random Numbers: Out of Uniform》