Matlabベクトル化二変量標準CDF

Question

ベクトル限界と相関rhoのためにcdfを計算する既存のコードがあるかどうか知っていますか？組み込みオプションmvncdfでは、行全体にわたってrhoを変更することはできません。

編集：

動作例では、各列xに対して同じミューおよびVを使用することである、すなわち

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

mvncdf(x,mu,sig); 

しかし、それぞれに異なるミュー及びSIGを入力する方法がありませんxの行。

Answer 1

これが必要な人は、完全にベクトル化されたバージョンがあります。非常に高い相関関係では精度は修正されていませんが、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 
%