加载中…
个人资料
supersasmacro
supersasmacro
  • 博客等级:
  • 博客积分:0
  • 博客访问:513,767
  • 关注人气:503
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
正文 字体大小:

SAS产生随机数的方法:随机数函数和CALL子程序

(2009-06-22 09:11:05)
标签:

sas

随机数

monte

carlo

蒙特卡罗

教育

分类: MonteCarlo蒙特卡罗

运用SAS进行Monte Carlo蒙特卡罗模拟(第四弹):

SAS产生随机数的方法:随机数函数和CALL子程序

 

本文未经作者同意严禁转载

 

 

1 随机数函数产生随机数序列

随机数函数产生随机数序列的语法:var=name(seed,<arg>),我们在前面的文章里使用的都是此类方法,变量var记录了由随机数种子为seed的随机函数name产生的一个随机数。我们举两个例子来深入说明随机数函数产生随机数序列的原理。

程序1:

DATA TEMP1(DROP=I);    

    DO I=1 TO 10      

       RUNI=RANUNI(123);

       SEED=RUNI*(2**31-1);      

       OUTPUT;    

    END;    

    RUN;

PROC PRINT DATA=TEMP1;

    RUN;

程序2:

DATA TEMP2(DROP=I);

    DO I=1 TO 5;

       RUNI1=RANUNI(123);

       RUNI2=RANUNI(456);***这里我们虽然指定了另一个随机数种子的值,但其实并不起作用;

       OUTPUT;

    END;

    RUN;

PROC PRINT DATA=TEMP2;

     RUN;

程序1的结果:

SAS产生随机数的方法:随机数函数和CALL子程序

 

程序2的结果:

SAS产生随机数的方法:随机数函数和CALL子程序

 

将程序1的结果和程序2的结果进行比较,程序2中的RUNI1是程序1中的RUNI的第1,3,5,7,9条数据,而程序2中的RUNI2是程序1中的RUNI的第2,4,6,8,10条数据。我们可以看到,虽然RUNI2=RANUNI(456);,但是这条语句中的随机数种子456并没有起作用。

RANUNI等随机函数的值产生后,SAS系统会隐含地将其转化为(0,1)的区间内,如果要看到这些随机数序列原始的值,可以通过程序1中的SEED=RUNI*(2**31-1);语句,结果如程序1输出结果中的seed所示。

 

2 CALL子程序产生随机数序列

CALL子程序产生随机数序列的语法如下:call name(seed,<arg>,var),与随机数函数产生随机数序列的语法类似。由上面的例子我们知道,随机数函数产生随机数序列时,其序列的值只由第一个随机数种子的值决定,而用CALL子程序时,每一次调用随机函数,都会重新产生新的随机数种子,下面举例进行讲解:

DATA TEMP3(DROP=I);

     RETAIN SEED1 123 SEED2 123 SEED3 123 SEED4 456;

     DO I=1 TO 10;

        RUNI1=RANUNI(SEED1);

        CALL RANUNI(SEED2,RUNI2);

        CALL RANUNI(SEED3,RUNI3);

        CALL RANUNI(SEED4,RUNI4);

        IF I=5 THEN DO; SEED1=456;***第五列开始,改变随机数种子的值;

                        SEED3=456;

                    END;

        OUTPUT;

     END;

     RUN;

PROC PRINT DATA=TEMP3;

     RUN;

结果如下:

SAS产生随机数的方法:随机数函数和CALL子程序

 

从上面的结果,我们可以看出:

前五行中,由于seed1=seed2=seed3,所以RUNI1,RUNI2和RUNI3的值相等,而SEED4=456,RUNI4与其它三个seed产生的值不一致,这就说明当使用CALL子程序进行随机数序列生成时,该随机数种子seed4起作用了,这是与使用函数生成随机数序列的一个重要的区别。

第六行开始,由于seed1和seed3的值有了变化,对随机数序列的变化情况如下:RUNI1和RUNI2的值并没有变化(可与前面两个程序的结果进行比对),RUNI3的值发行了变化。原因是在RUNI1=RANUNI(SEED1);这条语句中,SEED1的变化并不起作用,原因上一节已经讲了,因些RUNI1不发生变化,为什么RUNI3会变化呢,原因是CALL子程序调用函数生成随机数序列时,其随机数种子会更新,因此其随机数结果数据会发生变化。我们可以看到,RUNI3的第6到第10个随机数与RUNI4的第1到第5个随机数是一致的。

 

3 Seed为零时的随机数序列

前一篇文章已经讲过,随机数据集由随机数种子(seed)R0决定,它的值可以是(-(2^31-2),2^31-2)中的任何一个数。其中,当随机数种子(seed)的值大于0时,每次产生的R0的值都是一样的,但当随机数种子(seed)的值小于或等于零时,每次产生的R0的值是不同的,它会以系统的时间为种子产生随机数值。

当Seed为零时,每一次执行CALL子程序所产生的随机数序列结果都会不一致,下面举例进行讲解:

%MACRO SEED0;

 DATA TEMP4(DROP=I);

      SEED=0;

      DO I=1 TO 3;

         CALL RANUNI(SEED,RUNI);

         OUTPUT;

      END;

      RUN;

 PROC PRINT DATA=TEMP4;

      RUN;

%MEND;

TITLE 'First Macro Call';

%SEED0

TITLE 'Second Macro Call';

%SEED0

RUN;

结果如下:

First Macro Call

OBS SEED RANUI

1 21287539 0.00991

2 1972737807 0.91863

3 790022720 0.36788

Second Macro Call

OBS SEED RANUNI

1 156935322 0.07308

2 972755748 0.45297

3 2060025218 0.95927

我们可以看到,两个结果不一致。因此,我们可以将SEED=0看作是随机种子。

 

4 生成大量不重复的seed序列

在实际的应用中,我们经常会遇到需要大量随机数序列的情况,这时候我们就不能靠手工输入随要数种子。在第3节中,当SEED=0时,我们可以用这个随机种子产生大量的随机数序列,但是这里产生的随机数序列并不一定能保证这些随机数序列不重复。因此,本节介绍一个产生不重复的随机数种子的宏:

%MACRO SEEDGEN(FSEED=,LSTREAM=,NSEEDS=);

   DATA TEMP(KEEP=SEED);

        RETAIN SEED &FSEED;

        OUTPUT;

        DO I=1 TO MIN((&NSEEDS-1)*&LSTREAM,2**31-1) BY &LSTREAM;

           DO J=1 TO &LSTREAM;

              CALL RANUNI(SEED,X);

           END;

           OUTPUT;

        END;

        RUN;

   PROC PRINT DATA=TEMP;

        TITLE "List of &NSEEDS Seed Values";

        TITLE2 "Apart by &LSTREAM Numbers";

        RUN;

%MEND;

%SEEDGEN(FSEED=123,LSTREAM=20,NSEEDS=10);

其产生随机数种子的原理很简单,就是要产生多少个随机数种子,就按一定步长递增多少次,然后得到一个随机数作为种子。这个宏有个缺点,就是当步长*随机数种子数量>2**31-1时,可能得不到要求得到的随机数种子数量,原因是DO I=1 TO MIN (( &NSEEDS-1)* &LSTREAM, 2**31-1)  BY &LSTREAM;这条语句。但是这种实现方法比较简单,还是很有用处的。

 

参考资料

 

Xitao Fan, etc..Monte Carlo Studies: A Guide for Quantitative Researchers. SAS Institute Inc.,2002

 

本文章只用于学习,请不要用于商业目的,否则后果自负。

 

 

0

阅读 收藏 喜欢 打印举报/Report
  

新浪BLOG意见反馈留言板 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

新浪公司 版权所有