作者celestialgod (天)
看板MATLAB
标题Re: [讨论] 如何加快回圈产生数据
时间Tue Mar 8 10:18:38 2016
※ 引述《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://webptt.com/cn.aspx?n=bbs/MATLAB/M.1457403521.A.622.html
1F：推 popo14777: 谢谢大大 03/08 13:33
2F：→ popo14777: 不过小弟验证的时候，这两个还是有差耶，哪一个是正确 03/08 13:34
3F：→ popo14777: 的二元gamma呢? 03/08 13:34
4F：→ popo14777: 左边是大大你写的程式，右边是上网抓的 03/08 13:35
5F：→ popo14777: 还是就如你所说的没有正确的二元gamma? 03/08 13:36
都对，没有绝对定义的二元gamma
7F：→ sunev: 那第一个argument改1000就可以避掉回圈了 03/08 14:43
看起来是，我有GOOGLE到，不过没想到原PO他...
8F：推 sunev: @celestialgod，上两行是回popo14777的。我统计不好，没能 03/08 15:10
9F：→ sunev: 看出BivGamRND用的演算法，不过rho=0的话似乎只要用单变数 03/08 15:10
10F：→ sunev: 的gamma分布就好了？ 03/08 15:10
对，rho=0，用gmarnd就好，不过他应该是搭配他论文使用的
还是建议原PO自己去生成比较快，反正推导就这样而已
11F：推 popo14777: s大，第一个argument改1000? 听不太懂耶 03/08 15:24
12F：推 sunev: 试试看呗 03/08 16:15
13F：→ sunev: 第一个argument是指BivGamRND的第一个参数改成1000 03/08 17:17
14F：推 popo14777: 哦!s大 但我目的只是要在第一层回圈只产生一个数据而已 03/08 18:10
15F：→ popo14777: 因为我第一层回圈要算统计量所以只能丢一笔 03/08 18:10
16F：→ popo14777: 我试过在第一层与第二层之间产生1000组数据跟 03/08 18:11
17F：→ popo14777: 原先的方法第一层只跑产生一次,数据不太一样 03/08 18:13
18F：→ popo14777: 不知先1次产生1000笔在提出来好还是1次产生1笔好 03/08 18:14
19F：推 sunev: 看不懂你在说什麽，反正就你PO出来的东西而言，改1000最快 03/08 18:36
※ 编辑: celestialgod (140.109.74.87), 03/10/2016 12:13:34