手癢還是回一下好了
目標1:寫出U(0,1)
首先我們知道rand()是給出 0 ~ RAND_MAX
+--+--+...........................--+
0 1 2 RAND_MAX
全長為 RAND_MAX - 0 = RAND_MAX
所以產生U(0,1)的抽樣為
double Unf0(){
return (double)(rand())/(RAND_MAX));
}
Unf0在[0,1]內部點出的機率是對的
但在邊界上的機率會出問題(在0,1這兩點上)(consider open set)
在機率上0和1出現的機率為0 ie. p(0)=0 和 p(1)=0
為了避免這個問題可以把抽樣數加0.5
0.5 1.5 a b
| | | |
V V V V
+--+--+...........................--+----+
0 1 2 RAND_MAX RAND_MAX+1
a=RAND_MAX-0.5
b=RAND_MAX+0.5
全長為 RAND_MAX+1 - 0 = RAND_MAX + 1
double Unf1(){
return (double)(rand()+0.5)/(RAND_MAX+1.0));
}
但Unf1雖然沒有邊界機率為0的問題
但會有精度不夠的問題
這時候精度為 1 / (RAND_MAX+1)-->連float都不到
不過對於整數抽樣來說己經足夠了
其他的U(0,1)方法請參考 http://www.nrbook.com/a/bookcpdf.php 第7章
目標2:寫出U(a,b)
由[0,1]-->[a,b] 所以直接線性轉換就可以了
u~U(a,b)
v~U(0,1)
u=(b-a)*v+a;
程式碼如下
double Unfab(){
doubel rand=Unf1();
return (rand*(b-a)+a);
}
目標3:N個整數如何抽出來
1 2 3 4 5.....N
+--+--+--+--+-....+-
| | | | |
0.5 | 2.5 | ... |
1.5 3.5 (RAND_MAX+0.5)
原則:要讓每一個整數所量到的測度相同
也就是要抽k, 只要抽到[ k-0.5 ,k+0.5 ] 都算抽到 k
所以全空間應為[ 0.5 , N+0.5 ] 也就是 U(0.5 , N+0.5 )
以下例子對1,2,3,4,5個數做抽樣P次
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
int main(void){
srand(time(NULL));//#include<stdlib.h>, #include<time.h>
int P=100;
int a1[5]={0};
for(int i=0;i<P;i++){
//產生uniformly的分布
double a=(double)(rand()+0.5)/(RAND_MAX+1);//#include<stdlib.h>
double b=a*5+0.5;//轉成U(0.5,5.5)
//以下記錄每抽到多少
if(b<1.5){a1[0]++;}
else if(b<2.5){a1[1]++;}
else if(b<3.5){a1[2]++;}
else if(b<4.5){a1[3]++;}
else{a1[4]++;}
}
for(int i=0;i<5;i++){printf("%d\t",a1[i]);}//輸出
system("pause");//暫停功能
return 0;
}
[問題]為什麼不直接用%N就好了?
rand()%N : 產生 0 ~ N-1 的整數
rand() : 產生 0 ~ RAND_MAX 的整數
但 RAND_MAX / N 不一定是整除-->也就是有前幾個出現機率比較大
例子:
RAND_MAX=32767 對 5 個抽樣
RAND_MAX % 5 = 2;
也就是前3個抽到機率比較高
最後兩個機率比較低
所以才會去製造U(0,1)再轉回來對整數抽樣
這個問題可以參考 LPH66 的文章 文章代碼(AID): #1DSuIH6B
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.112.25.105
※ 文章網址: http://www.ptt.cc/bbs/C_and_CPP/M.1410237949.A.D40.html
※ 編輯: wope (140.112.25.105), 09/09/2014 12:51:12