用SAS模擬隨機(jī)數(shù)據(jù) 求PIE值

剛剛看到一本好書(shū)《統(tǒng)計(jì)模擬》作者叫羅斯[英文：Sheldon M. Ross. Simulation(4th Ed).Elsevier Inc..2006 ]. 顧名思義，這是一本描述怎么利用模擬一些符合統(tǒng)計(jì)學(xué)理論的數(shù)據(jù)，用途很廣，也就是說(shuō)實(shí)際中的任何數(shù)據(jù)的分布都符合某種統(tǒng)計(jì)學(xué)模型，于是在沒(méi)有得到真實(shí)數(shù)據(jù)之前，數(shù)據(jù)分析師可以通過(guò)模擬數(shù)據(jù)來(lái)研究這些現(xiàn)實(shí)中的問(wèn)題。如果通過(guò)模擬來(lái)研究未知問(wèn)題，可以說(shuō)得上是研究境界很高了?？偛荒苣玫揭恍?shí)際數(shù)據(jù)，畫(huà)個(gè)好看的圖，就覺(jué)得自己可畫(huà)遍天下了吧。

由于自己不是統(tǒng)計(jì)出生，但是受過(guò)統(tǒng)計(jì)學(xué)老師的循循教誨，凡事從簡(jiǎn)單開(kāi)始。于是goolge了一下，當(dāng)當(dāng)中有這本書(shū)的中文介紹：

本書(shū)系統(tǒng)闡述了統(tǒng)計(jì)模擬的一些實(shí)用方法和技術(shù)。在對(duì)概率的基本知識(shí)進(jìn)行了簡(jiǎn)單的回顧之后，介紹如何利用計(jì)算機(jī)產(chǎn)生隨機(jī)數(shù)以及如何利用這些隨機(jī)數(shù)產(chǎn)生任意分布的隨機(jī)變量、隨機(jī)過(guò)程等。然后討論了一些分析統(tǒng)計(jì)數(shù)據(jù)的方法和技術(shù)。如Bootstrap（自助法）、方差縮減技術(shù)等。接著講述了如何利用統(tǒng)計(jì)模擬來(lái)判斷所選的隨機(jī)模型是否擬合實(shí)際的數(shù)據(jù)。最后介紹MCMC及一些最新發(fā)展的統(tǒng)計(jì)模擬技術(shù)和論題，如隨機(jī)序列函數(shù)和隨機(jī)子集函數(shù)的評(píng)估。本書(shū)在每章的最后還提供了不同難度的習(xí)題。本書(shū)可作為高等院校數(shù)學(xué)、統(tǒng)計(jì)學(xué)、科學(xué)計(jì)算、保險(xiǎn)學(xué)、精算學(xué)等專業(yè)的教材，也可供工程技術(shù)人員和應(yīng)用工作者參考。

一看有很多不懂的術(shù)語(yǔ)，頓時(shí)心生敬仰，后面寫(xiě)著可以供“工程技術(shù)人員”參考，很顯然，我可以是這本書(shū)的讀者。全書(shū)近300頁(yè)，在今天剩下不多的時(shí)間內(nèi)，很顯然，我決定不去看這本書(shū)，哪怕是一個(gè)字。還是老師的教導(dǎo)，從簡(jiǎn)單開(kāi)始。google到一個(gè)有趣的問(wèn)題：用統(tǒng)計(jì)模擬計(jì)算圓周率pie值。Forcode提供了一種用excel求解pie的詳細(xì)過(guò)程，然后有人用Mathematica計(jì)算出來(lái)了。Hujiangtang很仔細(xì)的闡述了什么是隨機(jī)數(shù)？什么是蒙特卡羅模擬？為什么選擇SAS做蒙特卡羅模擬？SAS怎么做出來(lái)這些，還可以做哪些哪些分布…… 其中引用這個(gè)用統(tǒng)計(jì)模擬計(jì)算pie的例子和上述用excel和mathematica的方法，可是就是不提供解決這個(gè)問(wèn)題的SAS代碼，在我看來(lái)，這是不可想象，我很仔細(xì)的找了半天，未果，于是我自己編了下列代碼，用SAS來(lái)實(shí)現(xiàn)模擬pie值。

借用別人寫(xiě)的原理和圖來(lái)說(shuō)明一下，

原理很簡(jiǎn)單：

1）生成隨機(jī)數(shù)——生成n個(gè)均勻落在正方形內(nèi)的點(diǎn)；

2）對(duì)落在正方形內(nèi)的n個(gè)點(diǎn)，數(shù)一數(shù)正好落在圓里面的點(diǎn)的個(gè)數(shù)，假設(shè)為k（另外n-k個(gè)點(diǎn)就落在圓外面的正方形區(qū)域內(nèi)）。數(shù)據(jù)分析師培訓(xùn)

3）k/n就可以大致認(rèn)為是圓的面積與正方形的面積之比，另其等于pai/4，就可以求出圓周率∏的估計(jì)值。n越大，算出來(lái)的pai-值越精確

SAS實(shí)現(xiàn)代碼：

%let num_s=10000;
data ex;
do i= 1 to &num_s.;
x=ranuni(0); y=ranuni(0);
area=sqrt(x**2+y**2);
output;
end;
run;

proc sort;
by area;
run;

data ex2;
set ex;
y1=.;y2=.;
num=_n_;
if area<=1 then y1=y ; else y2=y;
if area>1 and lag(area)<=1 then
do;
pie=(num*4)/&num_s.;
sas_pie=constant(‘pi’);
put pie= sas_pie=;
end;
run;
FILENAME file “c:\simulation.png”;
goptions reset=all hsize=8cm vsize=8cm  noborder device=png gsfname=file;;
symbol1 v=dot color=’red’ height=0.3;
symbol2 v=dot color=’blue’ height=0.3;
proc gplot;
plot  y1*x y2*x/overlay noaxis;
run;
quit;

模擬的點(diǎn)數(shù)為 100時(shí)pie=3；1000時(shí)，pie=3.068； 10000時(shí)， pie=3.1392；100000時(shí)pie=3.13543；1000000時(shí)，pie=3.141524，10000000時(shí)，pie=3.1418008。當(dāng)然SAS也提供了pie的精確值供使用，使用函數(shù)constant就可以：sas_pie=constant(‘pi’); 值為：3.1415926536。上面的示例只是為了演示統(tǒng)計(jì)模擬的使用，其實(shí)SAS提供了很多種隨機(jī)函數(shù)，幾乎所有分布的隨機(jī)函數(shù)供大家使用，這里就不詳述了.CDA數(shù)據(jù)分析師培訓(xùn)