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)]