推 romotting:謝謝~我再研究一下(我不太會寫程式) 迭代那邊是我搞錯了 11/23 17:38
→ romotting:抱歉> < 11/23 17:38
推 DEATHX:E大出品，必屬精品。 11/23 21:24
查了一下資料，其實真的不難，由於 internet 掛了，所以這篇會放 code
( 沒辦法用置底文放 code)，等好的時候再補上置底文連結。
※ 引述《romotting ()》之銘言：
: 開發平台(Platform):DevC++
: 額外使用到的函數庫(Library Used): (Ex: OpenGL, ...)
: 問題(Question)：
: 1.原本我想要讓n1(AA基因型):n2(Aa基因型):n3(aa基因型)個數隨機服從multinomial
: distrbution,但是我只google到均勻分配的,所以只好假設n=6000,p=0.6 >>n1=2160
: n2=2880,n3=960,請問有神可以告訴我怎麼隨機生成多項分配嗎??
http://en.wikipedia.org/wiki/Multinomial_distribution
wiki 上拉到 To simulate a multinomial distribution，
它有範例跟你講，用 excel 示意，我沒意會錯的話它的意思如下 ( 拿你的例子講 )
(1) 先決定各基因機率值 P(n1) = 0.36, P(n2) = 0.48, P(n3) = 0.16
初值化 n1=n2=n3=0 , 令母體數 = N, 最後 n1+n2+n3=N
(2) 據 P(n1~n3), 計算累計機率 CP(n1)=0.36 , CP(n2)=0.84, CP(n3)=1.0
(3) 接下來是隨機式輪盤法時間。
for(i=0; i<N; ++i) {
(3.1) 隨機產生一個 [0, 1) 亂數 rnd ,
(3.2) 找 j, 滿足 CP(nj) <= rnd < CP(nj+1)
}
如果你有 K 過基因演算法的話，就是隨機式輪盤法(選擇運算子)
範例碼如下
const size_t PCnt = 6000; /* 母體個數 */
const size_t N = 3;
const double Prob[N] = {0.36, 0.48, 0.16};
double ColumnProb[N] ;
size_t i, CCnt[N] = {0}; /* CCnt[i] : 子代 i 個數 */
size_t Interval;
double rnd;
// 計算累計機率
ColumnProb[0] = Prob[0];
for(i=1; i<N; ++i)
ColumnProb[i] = Prob[i] + ColumnProb[i-1];
// 進行分佈
for(i=0; i<PCnt; ++i) {
rnd = rand() / (RAND_MAX + 1.0);
// Find Interval
for(Interval=0; Interval<N; ++Interval)
if(ColumnProb[Interval] >= rnd)
break;
++CCnt[Interval];
}
最後 CCnt 就是你要的。
: 2.我產生出來的隨機機率p想要固定住,所以寫p0==p;之後想看遞迴後和p0絕對值差10的
: -8次,但是好像有問題,因為結果出來一堆0,我在想是不是因為p0==p;那如果想要讓產
: 生出來的亂數值固定要怎麼設?
<恕刪>
: p=rand()/(RAND_MAX + 1.0);
: p0==p;
: while(i<1000)
: {
: p=((n1+n2)/6000)*(1/2-p);
: if(abs(p0-p)<=10^-8)
: cout << p << endl;
: i++;
: }
: system("pause");
: return 0;
: }
這裡是我猜的。
我猜你對於「收斂」的定義有問題。一般而言，類似的演算法大概是，
先隨便找一個點 (p0) 當初點，然後不停迭代，第一次迭代結果得到 p1；
第二次迭代結果得到 p2，...，第 n 次迭代結果得到 pn，
收斂應該指的是 abs( pn - pn+1) <= EPS 視為收斂，
而非 abs(pn - p0) <= EPS，所以 我猜 大概是要改這樣
double p = rand() / (RAND_MAX+1.0);
double p0 = p;
const double EPS = 1E-8; // 收斂誤差
const size_t MAX_ITER = 1000U; // 最大迭代次數
for(i=0; i<MAX_ITER; ++i){
p = (CCnt[0]+CCnt[1]) * (0.5 - p) / PCnt;
if( abs(p-p0) < EPS) // 已達終止條件
break;
p0 = p;
}
cout << " p = " << p << endl;
cout << "iterator = " << i << endl;
和你寫的有些不同，但注意的是，
一些 magic number 請養成習慣寫成 constant variable,
以後寫書面、維護程式時絕對有幫助。
上述 code 我測的確會收斂，普遍是在 100 多次，收斂到 0.22~0.25 不等。
< 以下吃光光 >
--
就算把新鮮的肝拿回去，還是一樣寫碼到禿頭，加班到天亮，
永遠當老闆的傀儡 你是不是想這麼做？
是的話你就拿回去~ 拿啊！！
九世宅男 : 下輩子不要再讓我幹工程師了 ~
< Kuso 星爺語錄 >
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 180.177.76.161