Re: [討論] 如何加快迴圈產生數據

[celestialgod (天)]

※ 引述《popo14777 (草草)》之銘言：
: 小弟想要跑三層迴圈的ARL，以下是我的程式碼
: tic
: for i=1:1000
:   mvnrnd([0 0],[1 0; 0 1]);
: end
: toc
: tic
: for i=1:1000
:   [x1 x2]=BivGamRND(1, 4, 1, 4, 1, 0);
:     Q=[x1 x2];
: end
: toc
: 結果為
: Elapsed time is 0.056313 seconds.
: Elapsed time is 50.921110 seconds.
: 第一種產生二元常態與第二種產生二元Gamma差了1000倍左右...
: 這只是第一層而已，第二層j要重複1000次，第三層k要run 100次....
: 想請問大大如何讓我的二元gamma產生數據快一點呢?
: 謝謝

因為雙元的GAMMA沒有統一的規定 (其實像是exponential, poisson都是)

如果你覺得他的生的慢，你可以考慮自己寫

利用線性組合的方式，只是必須下一點功夫去算covariance


X1 ~ Gamma(a1, b1),  X2 ~ Gamma(a2, b2),  X1 is independent of X2

Let Y1 = X1 + p * X2, Y2 = p * X1 + X2

E(Y1) = a1 * b1 + p * a2 * b2

E(Y2) = p*a1*b1 + a2*b2

Var(Y1) = Var(X1) + p^2 * Var(X2) = a1*b1^2+a2*(b2*p)^2

Var(Y2) = a1 * (p*b1)^2 + a2*b2^2

Cov(Y1, Y2) = Cov(X1 + p*X2, p*X1+X2) = p*Var(X1) + p*Var(X2) + 2*p*Cov(X1,X2)

= p*a1*b1^2 + p*a2*b2^2  (Cov(X1,X2) = 0 because X1 is indep. of X2)

推完之後，你就可以控制你的p，去控制共變異數的大小 (或是相關係數)

然後生成的時候 就是


function Q = bivgmarnd(N, a1, b1, a2, b2, p)

X1 = gamrnd(a1,b1, N, 1);
X2 = gamrnd(a2,b2, N, 1);

Q = [X1, X2] * [1, p; p, 1];

這樣生就會快很多


想要控制rho可以這樣寫：

好讀版：http://pastebin.com/wTP3qTRJ

function Q = bivgmarnd(N, a1, b1, a2, b2, rho)
syms p;
p_ = double(solve((p*a1*b1^2 + p*a2*b2^2) / (sqrt(a1*b1^2+a2*(b2*p)^2) *
sqrt(a1 * (p*b1)^2 + a2*b2^2)) == rho, p));
rho_ = (p_*a1*b1^2 + p_*a2*b2^2) ./ (sqrt(a1*b1^2+a2*(b2*p_).^2) .* sqrt(a1 *
(p_*b1).^2 + a2*b2^2));
validate = find(abs(rho_ - rho) < 1e-3);
if isempty(validate)
  error('cannot find the real solution for (p1, p2)')
else
  p_ = p_(validate(1));
end

X1 = gamrnd(a1, b1, N, 1);
X2 = gamrnd(a2, b2, N, 1);

Q = [X1, X2] * [1, p_; p_, 1];
end

corr(bivgmarnd(1e6, 3, 2, 3, 1, 0.6))
% ans =
%     1.0000    0.6015
%     0.6015    1.0000
corr(bivgmarnd(1e6, 3, 2, 3, 1, 0.9))
% ans =
%     1.0000    0.9001
%     0.9001    1.0000

--
R資料整理套件系列文：
magrittr    #1LhSWhpH (R_Language) http://tinyurl.com/1LhSWhpH
data.table  #1LhW7Tvj (R_Language) http://tinyurl.com/1LhW7Tvj
dplyr(上)   #1LhpJCfB (R_Language) http://tinyurl.com/1LhpJCfB
dplyr(下)   #1Lhw8b-s (R_Language)
tidyr       #1Liqls1R (R_Language) http://tinyurl.com/1Liqls1R

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.109.74.87
※ 文章網址: https://www.ptt.cc/bbs/MATLAB/M.1457403521.A.622.html
推 popo14777: 謝謝大大
→ popo14777: 不過小弟驗證的時候，這兩個還是有差耶，哪一個是正確
→ popo14777: 的二元gamma呢?
→ popo14777: 左邊是大大你寫的程式，右邊是上網抓的
→ popo14777: 還是就如你所說的沒有正確的二元gamma?
都對，沒有絕對定義的二元gamma

推 sunev: 你原來的 BivGamRND 是從這來的嗎？http://0rz.tw/koIZB
→ sunev: 那第一個argument改1000就可以避掉迴圈了
看起來是，我有GOOGLE到，不過沒想到原PO他...

推 sunev: @celestialgod，上兩行是回popo14777的。我統計不好，沒能
→ sunev: 看出BivGamRND用的演算法，不過rho=0的話似乎只要用單變數
→ sunev: 的gamma分布就好了？
對，rho=0，用gmarnd就好，不過他應該是搭配他論文使用的

還是建議原PO自己去生成比較快，反正推導就這樣而已
推 popo14777: s大，第一個argument改1000? 聽不太懂耶
推 sunev: 試試看唄
→ sunev: 第一個argument是指BivGamRND的第一個參數改成1000
推 popo14777: 哦!s大 但我目的只是要在第一層迴圈只產生一個數據而已
→ popo14777: 因為我第一層迴圈要算統計量所以只能丟一筆
→ popo14777: 我試過在第一層與第二層之間產生1000組數據跟
→ popo14777: 原先的方法第一層只跑產生一次,數據不太一樣
→ popo14777: 不知先1次產生1000筆在提出來好還是1次產生1筆好
推 sunev: 看不懂你在說什麼，反正就你PO出來的東西而言，改1000最快
※ 編輯: celestialgod (140.109.74.87), 03/10/2016 12:13:34