Matlab vectorized bivariate標準正常CDF

有誰知道是否有現有的代碼來計算cdf的向量限制和相關性rho？內置選項mvncdf不允許rho在行之間變化。Matlab vectorized bivariate標準正常CDF

編輯：

的作品是每行的X相同畝和V的一個例子，即

x = rand(10,2); 
mu = [0 0]; 
sig = [1 .5; .5 1]; 

mvncdf(x,mu,sig);

但是，有沒有辦法進入一個不同畝，SIG每個x的行。

來源

2017-06-18 user2276896

@SardarUsama Matlab循環是最後的手段。 – user2276896

您可以使用Alan Genz提供的[m-file]（http://www.math.wsu.edu/faculty/genz/software/matlab/tvtl.m）。函數'bvnu'可以有效地矢量化。所以你應該在使用之前用矢量化的形式重寫它。有關更多參考資料，請閱讀本文[論文]（http://www.math.wsu.edu/faculty/genz/papers/bvnt.pdf） – rahnema1

如果有人需要這個，這是一個完全向量化的版本。對於非常高的相關性，它的準確性沒有修改，但對於rho < .9應該沒有問題。

function p = bvnl(dh, dk, r) 
%BVNL 
% A function for computing bivariate normal probabilities. 
% bvnl calculates the probability that x < dh and y < dk. 
% parameters 
%  dh 1st upper integration limit 
%  dk 2nd upper integration limit 
%  r correlation coefficient 
% Example: 
% p = bvnl(3, 1, .35) 
% 

% Original author Alan Genz, modified to be vectorized. 

     p = bvnu(-dh, -dk, r); 
% 
% end bvnl 
% 
function p = bvnu(dh, dk, r) 
%BVNU 
% A function for computing bivariate normal probabilities. 
% bvnu calculates the probability that x > dh and y > dk. 
% parameters 
%  dh 1st lower integration limit 
%  dk 2nd lower integration limit 
%  r correlation coefficient 
% Example: p = bvnu(-3, -1, .35) 
% Note: to compute the probability that x < dh and y < dk, 
%  use bvnu(-dh, -dk, r). 
% 

    if dh == inf | dk == inf, 
     p = 0; 
    elseif dh == -inf, 
     if dk == -inf, 
      p = 1; 
     else 
      p = phid(-dk); 
     end 
    elseif dk == -inf, 
     p = phid(-dh); 
    elseif r == 0, 
     p = phid(-dh)*phid(-dk); 
    else 
     tp = 2*pi; 
     h = dh; 
     k = dk; 
     hk = h.*k; 
     bvn = 0; 
    if abs(r) < 0.3  % Gauss Legendre points and weights, n = 6 
     w(1:3) = [0.1713244923791705 0.3607615730481384 0.4679139345726904]; 
     x(1:3) = [0.9324695142031522 0.6612093864662647 0.2386191860831970]; 
    elseif abs(r) < 0.75 % Gauss Legendre points and weights, n = 12 
     w(1:3) = [.04717533638651177 0.1069393259953183 0.1600783285433464]; 
     w(4:6) = [0.2031674267230659 0.2334925365383547 0.2491470458134029]; 
     x(1:3) = [0.9815606342467191 0.9041172563704750 0.7699026741943050]; 
     x(4:6) = [0.5873179542866171 0.3678314989981802 0.1252334085114692]; 
    else     % Gauss Legendre points and weights, n = 20 
     w(1:3) = [.01761400713915212 .04060142980038694 .06267204833410906]; 
     w(4:6) = [.08327674157670475 0.1019301198172404 0.1181945319615184]; 
     w(7:9) = [0.1316886384491766 0.1420961093183821 0.1491729864726037]; 
     w(10) = 0.1527533871307259; 
     x(1:3) = [0.9931285991850949 0.9639719272779138 0.9122344282513259]; 
     x(4:6) = [0.8391169718222188 0.7463319064601508 0.6360536807265150]; 
     x(7:9) = [0.5108670019508271 0.3737060887154196 0.2277858511416451]; 
     x(10) = 0.07652652113349733; 
    end, w = [w w]; x = [1-x 1+x]; 

    hs = (h.*h + k.*k)/2; 
    asr = asin(r)/2; 
    sn = sin(asr*x); 
    bvn = exp(bsxfun(@minus,bsxfun(@times,sn,hk),hs)./(1-sn.^2))*w'; 
    bvn = bvn.*asr./tp + phid(-h).*phid(-k); 

    p = max(0, min(1, bvn)); 
    end 
% 
% end bvnu 
% 
function p = phid(z), p = erfc(-z/sqrt(2))/2; % Normal cdf 
% 
% end phid 
%

來源

2017-06-19 16:42:58 user2276896

Matlab vectorized bivariate標準正常CDF

回答

相關問題