Jiangtang's profile技止于此BlogListsNetwork Tools Help

Blog


    1/22/2008

    SAS和蒙特卡罗模拟(3):SAS随机数函数及CALL子程序

    **************************************************************************************************************************

    SAS和蒙特卡罗模拟(1):开篇》——简介,通过例子建立起蒙卡的直观概念;参考软件包及书目

    SAS和蒙特卡罗模拟(2):随机数基础》——随机数入门;参考书目

    **************************************************************************************************************************************************

    一、SAS随机数函数和CALL子程序 

    SAS系统产生随机数,两种方式,利用SAS函数(Functions)或者CALL子程序(CALL Routines),它们的语法格式是(具体的区别容后讨论):

              方式           代码                         说明
    函数 var=name(seed,<arg>) var为存储随机数列的变量,name为特定的分布函数形式,seed为随机数种子,<arg>为特定分布要求的参数(可选)
    CALL子程序 call name(seed,<arg>,var) 同上,记得seed=0, ±1,±2, , ± (2**31-2)

    SAS可用的随机数函数列表如下(可以参见SAS Help and Documentation-SAS Products-Base SAS-SAS Language Dictionary-Functions and CALL Routines-Functions and CALL Routines by Category):

            分布     SAS函数         说明
    二项分布(Binomial) ranBin(seed,n,p) n为独立实验的次数,p为成功概率
    指数分布(Exponential) ranExp(seed)  
    正态分布(Normal) ranNor(seed),or normal(seed) ranNor和normal是同质的,但normal没有相对应的CALL子程序
    泊松分布(Poisson) ranPoi(seed,m) m为均值
    均匀分布(Uniform) ranUni(seed),or uniform(seed) ranUni和uniform是同质的,但uniform没有相对应的CALL子程序
    柯西分布(Cauchy) ranCau(seed)  
    伽玛分布(Gamma) ranGam(seed,a) a>0为形状参数
    由分布律表格决定的离散分布(tabled probability distribution) ranTbl(seed,p1,p2,...pn) ∑p=1
    三角分布(Triangular) ranTri(seed,h) h为斜边(最可能值)

    上面的随机函数,除了normal和uniform,都可以由直接相应的CALL子程序调用。

    二、SAS随机数函数和CALL子程序:比较

    用SAS随机数函数同时创建的多个随机数变量其实都属于同一个随机数列。这话费解,一个例子先,创建两个随机数变量,各包含3个记录,其中x1的种子为123,x2的种子为456:

    data ranuni(drop=i);
        retain seed1 123 seed2 456;
        do i=1 to 3;
            x1=ranuni(seed1);
            x2=ranuni(seed2);       
            output;
        end;
    run;
    proc print data=ranuni;run;

    结果为:

    Obs    seed1    seed2       x1         x2

    1      123      456     0.75040    0.32091
    2      123      456     0.17839    0.90603
    3      123      456     0.35712    0.22111

    好像没什么异样。我们把上面的x1增加为6个记录:

    data ranuni2(drop=i);
        retain seed1 123;
        do i=1 to 6;
            x1=ranuni(seed1);
            output;
        end;
    run;
    proc print data=ranuni2;run;

    结果如下,把上下用红色标出的数字对照看一看:

    Obs    seed1       x1

    1      123     0.75040
    2      123     0.32091
    3      123     0.17839
    4      123     0.90603
    5      123     0.35712
    6      123    0.22111

    什么意思?在第一段代码中,x2的种子根本不起作用,把x2的记录安插到x1里,看起来就是用种子123产生的随机数列加长了而已。x2的第一个值并不是由种子456产生的,而是产生第一个x1后的下一个值,x1、x2其实属于同一个随机数列,尽管x2的种子被指定为456,而x1的被指定为123。现在就可以重复上面的一句话:用SAS随机数函数同时创建的多个随机数变量其实都属于同一个随机数列

    用CALL子程序就能够同时产生独立的随机数列。

    data ranuni3(drop=i);
        retain seed3 123 seed4 456;
        do i=1 to 3;
            call ranuni(seed3,x3);
            call ranuni(seed4,x4);
            output;
        end;
    run;
    proc print data=ranuni3;run;

    结果如下:

    Obs       seed3        seed4         x3         x4

    1     1611463328    736440516   0.75040    0.34293
    2      689153326    774069794    0.32091    0.36045
    3      383088854    686944750    0.17839    0.31988

    以上x3就是x1。x1和x3的初始种子都是123,但x3那个结果显示的种子是当前种子值。要在SAS随机数函数语句中显示随机种子的当前值,可以由以下代码给出:

    data ranuni4(drop=i);
        retain seed1 123;
        do i=1 to 3;
            x1=ranuni(seed1);
            seed=x1*(2**31-1);   
            output;
        end;
    run;
    proc print data=ranuni4;
        var seed1 seed x1;
    run;

    结果如下,可以跟上面由CALL子程序得出的结果对照:

    Obs    seed1       seed          x1

    1      123     1611463328   0.75040
    2      123      689153326    0.32091
    3      123      383088854    0.17839

    ---------参考资料---------

    1. Xitao Fan, etc..Monte Carlo Studies: A Guide for Quantitative Researchers. SAS Institute Inc.,2002
    2. 朱世武《SAS编程技术与金融数据处理》,北京:清华大学出版社,2003

    Comments (2)

    Please wait...
    Sorry, the comment you entered is too long. Please shorten it.
    You didn't enter anything. Please try again.
    Sorry, we can't add your comment right now. Please try again later.
    To add a comment, you need permission from your parent. Ask for permission
    Your parent has turned off comments.
    Sorry, we can't delete your comment right now. Please try again later.
    You've exceeded the maximum number of comments that can be left in one day. Please try again in 24 hours.
    Your account has had the ability to leave comments disabled because our systems indicate that you may be spamming other users. If you believe that your account has been disabled in error please contact Windows Live support.
    Complete the security check below to finish leaving your comment.
    The characters you type in the security check must match the characters in the picture or audio.

    To add a comment, sign in with your Windows Live ID (if you use Hotmail, Messenger, or Xbox LIVE, you have a Windows Live ID). Sign in


    Don't have a Windows Live ID? Sign up

    No namewrote:
    您需要二手液晶显示屏废旧液晶屏么?我们是不折不扣的二手液晶屏、旧液晶屏大批发商,长期大量供应可再利用的旧液晶屏。我公司提供的各种尺寸的二手液晶屏, 不同厚薄如笔记本屏,均已经过我们严格的分类,检验,测试流程。请访问协力液晶屏www.sceondhandlcd.com[ghfdeeaeaecfgca]
    Nov. 21
    Nov. 9

    Trackbacks

    The trackback URL for this entry is:
    http://johnthu.spaces.live.com/blog/cns!2053CD511E6D5B1E!379.trak
    Weblogs that reference this entry
    • None