Re: Inconsistency in different ways of calculating CDF of bivariate

• To: mathgroup at smc.vnet.net
• Subject: [mg106433] Re: Inconsistency in different ways of calculating CDF of bivariate
• From: Ray Koopman <koopman at sfu.ca>
• Date: Tue, 12 Jan 2010 04:50:26 -0500 (EST)
• References: <higdkt\$kgl\$1@smc.vnet.net>

```On Jan 11, 3:53 pm, Stephen <aultman.step... at gmail.com> wrote:
> I would like to calculate the CDF of a bivariate normal in some C++
> code. My first thought was to reduce the dimensionality of the
> integral by using the conditonal distribution of the bivariate normal
> when one variable is known (see link below). When I try to calculate
> the CDF  this way I get an incorrect result. In the code pasted below
> I defined a function, temp1, that uses mathematica's built in
> bivaraite normal CDF function. The function temp2 attempts to
> calculate the same quantity using a one dimension integral and the
> conditional distribution of the bivariate normal.
>
> It is more than just rounding errors. For some values of rho, the
> error is as much as 0.03, which is a lot in probability terms. Mostly
> I am curious why this doesn't work from a mathematical statistics
> perspective. There are ways to code around this, but I am curious if
> there is a bigger math issue at play. Any thoughts?
>
> mu1 = 0; mu2 = 0; sigma1 = 1; sigma2 = 1; rho = .1; z1 = 0; z2 = 0;
> Needs["MultivariateStatistics`"]
> temp1[mu1_, mu2_, sigma1_, sigma2_, rho_, z1_, z2_] :=
>  CDF[MultinormalDistribution[{mu1,
>     mu2}, {{sigma1^2, rho*sigma1*sigma2}, {rho*sigma1*sigma2,
>      sigma2^2}}], {z1, z2}]
> temp2[mu1_, mu2_, sigma1_, sigma2_, rho_, z1_, z2_] :=
>  NIntegrate[
>   PDF[NormalDistribution[mu1, sigma1], x]*
>    CDF[NormalDistribution[
>      mu2 + (sigma2/sigma1)*rho (x - mu1), (1 - rho^2)*sigma2^2],
>     z2], {x, -Infinity, z1}]
>
> Plot[temp1[mu1, mu2, sigma1, sigma2, rho, z1, z2] -
>   temp2[mu1, mu2, sigma1, sigma2, rho, z1, z2], {rho, -.99, .99}]
>
> http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Conditional_distributions

In temp2, change  (1 - rho^2)*sigma2^2  to  Sqrt[1 - rho^2]*sigma2 .
(You want the standard deviation, not the variance.)

Here are two other routines that you might want to try. They use
equations from:  Pearson, K. (1901). Mathematical contributions to
the theory of evolution - VII: On the correlation of characters not
quantitatively measurable. Philosophical Transactions of the Royal
Society of London, Series A, 195, 1-47.

x = (z1 - mu1)/sigma1, y = (z2 - mu2)/sigma2, r = rho.
cdf get the lower left area, < x and < y; ccdf gets the upper right
area, > x and > y (Abramowitz & Stegun's "L").

cdf[x_,y_,r_] := .25 Erfc[-Sqrt[.5]*x] Erfc[-Sqrt[.5]*y] +
Block[{u = -.5(x^2 + y^2), v = N[x*y]}, Which[u == 0, ArcSin[r],
v == 0, NIntegrate[Exp[u/Cos[t]^2],{t,0,ArcSin@r}],
True, NIntegrate[Exp[(u + v*Sin[t])/Cos[t]^2],{t,0,ArcSin@r}]]/(2 Pi)]

ccdf[x_,y_,r_] := .25 Erfc[Sqrt[.5]*x] Erfc[Sqrt[.5]*y] +
Block[{u = -.5(x^2 + y^2), v = N[x*y]}, Which[u == 0, ArcSin[r],
v == 0, NIntegrate[Exp[u/Cos[t]^2],{t,0,ArcSin@r}],
True, NIntegrate[Exp[(u + v*Sin[t])/Cos[t]^2],{t,0,ArcSin@r}]]/(2 Pi)]

```

• Prev by Date: Re: restricting interpolating functions to be positive
• Next by Date: Re: Re: Re: More /.{I->-1} craziness
• Previous by thread: Inconsistency in different ways of calculating CDF of bivariate
• Next by thread: syntax extension