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

看板MATLAB作者celestialgod (天)時間8年前 (2016/03/08 10:18)推噓7(7推 0噓 12→)

留言19則, 2人參與討論串2/2 (看更多)

※ 引述《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

03/08 13:33, , 1^F

03/08 13:33, 1^F

→

popo14777

03/08 13:34, , 2^F

03/08 13:34, 2^F

→

popo14777

03/08 13:34, , 3^F

03/08 13:34, 3^F