推 VictorTom:太強了~~拜....<(_ _)> 05/05 00:49
※ 引述《inanity (陰晴圓缺)》之銘言：
: 開發平台(Platform): (Ex: VC++, GCC, Linux, ...)
: C++
: 額外使用到的函數庫(Library Used): (Ex: OpenGL, ...)
: 問題(Question)：
: 看不懂網路上的程式碼，尤其是在while迴圈內外都有一個kappa變數，不知道該怎樣解釋
: 餵入的資料(Input)：
看起來 input 有這些 (假設 type 都是 double)
scalar(double): mu_s, mu_m, hlp, k_c, th_c, sv_c, ro_c, n_sim
array(double[]): v, j, ky, kv, y
pointer(double*): v_pre
在 c 裡 pointer 可以指到 array 的某個元素
例如 v_it 是一個 pointer，v_it = v 表示 v_it 指到 v 這個 array 的開頭
然後因為 c 的 index 是從 0 開始算的，所以 *v_it 就是 v[0]
之後 v_it++ 表示 v_it 現在指向下一個了 (所以是 v[1])
所以要改成 matlab 主要就是下面幾點：
* 把 c 的 array 改成 matlab 的 vector，pointer 改成 int
c version matlab version
v_it = v v_it = 1
*v_it v(v_it)
v_it++ v_it = v_it + 1
* matlab 的 index 是從 1 開始算，跟 c 從 0 開始算會差一個
* v_pre 在 c 裡應該是指向 v 的某個元素
在 matlab 裡要先算出是那一個
* pow(a, b) -> a^b
* gsl_ran_gaussian，在 matlab 裡要找到相同功能的函式取代
(下面的 code 我是用 randn，要確認一下是不是完全相同)
: 預期的正確結果(Expected Output)：
: 錯誤結果(Wrong Output)：
: 程式碼(Code)：(請善用置底文網頁, 記得排版)
: #include<math.h>
: #include "bayes.h"
: kappa1 = 1.0/pow(mu_s,2.0);
: kappa2 = mu_m/pow(mu_s,2.0);
: v_it=v; j_it=j; ky_it=ky; kv_it=kv; y_it=y;
: while(v_it<v_pre) {
: hlp = *v_it++;
: e_v = ((*v_it) - hlp - k_c * (th_c - hlp) - (*kv_it)*(*j_it))...
: /sqrt(sv_c);
: e_y = (*y_it) - (*ky_it)*(*j_it) - ro_c*e_v;
: kappa2 += e_y/(1.0 - pow(ro_c,2.0))/hlp;
: kappa1 += 1/hlp/(1.0 - pow(ro_c,2.0));
: ++j_it; ++ky_it; ++kv_it; ++y_it;
: }
: hlp = *v_it;
: e_y = (*y_it) - (*ky_it)*(*j_it);
: kappa2+=e_y/hlp;
: kappa1+=1/hlp;
: kappa1=1/kappa1;
: mu_c = gsl_ran_gaussian(rv_g,sqrt(kappa1))+ kappa1 * kappa2;
: Emu += mu_c/double(nsim);
kappa1 = 1.0/(mu_s^2.0);
kappa2 = mu_m/(mu_s^2.0);
v_it = 1; j_it = 1; ky_it = 1; kv_it = 1; y_it = 1;
while v_it < v_pre
hlp = v(v_it);
v_it = v_it + 1;
e_v = v(v_it) - hlp - k_c*(th_c-hlp) - kv(kv_it)*j(j_it)/sqrt(sv_c);
e_y = y(y_it) - ky(ky_it)*j(j_it) - ro_c*e_v;
kappa2 = kappa2 + e_y/(1.0 - (ro_c^2.0))/hlp;
kappa1 = kappa1 + 1/hlp/(1.0 - (ro_c^2.0));
j_it = j_it + 1; ky_it = ky_it + 1; kv_it = kv_it + 1; y_it = y_it + 1;
end
hlp = v(v_it);
e_y = y(y_it) - ky(ky_it)*j(j_it);
kappa2 = kappa2 + e_y/hlp;
kappa1 = kappa1 + 1/hlp;
kappa1 = 1/kappa1;
mu_c = sqrt(kappa1)*randn(1) + kappa1*kappa2;
Emu = Emu + mu_c/double(nsim);
我照著翻而已，沒測過
加減參考吧...
: 補充說明(Supplement)：
: 小弟我學的是MATLAB，所以對C++不熟，想要把這個程式碼改成MATLAB語言
: 如果板友對這個程式碼還需要其他補充資料，請告知
: 謝謝各位大大
: PS 這是ERAKER 2003年MCMC的程式碼
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.112.30.49