Re: [R] distribution of the product of two correlated normal
Yu, Xuesong schrieb:
Many thanks to Peter for your quick and detailed response to my question.
I tried to run your codes, but seems like u is not defined for functions fp
and fm. what is u?
I believe t=X1*X2
nen0 - m2+c0*u ## for all u's used in integrate: never positive
no, this is not the problem; u is the local integration variable
in local functions f, fm, fp over which integrate() performs
integration;
it is rather the eps = eps default value passed in functions
f, fm, fp which causes a recursive default value reference - problem;
change it as follows:
###
#code by P. Ruckdeschel, [EMAIL PROTECTED], rev. 04-25-06
###
#
#pdf of X1X2, X1~N(m1,s1^2), X2~N(m2,s2^2), corr(X1,X2)=rho, evaluated at t
#
# eps is a very small number to catch errors in division by 0
###
#
dnnorm - function(t, m1, m2, s1, s2, rho, eps = .Machine$double.eps ^ 0.5){
a - s1*sqrt(1-rho^2)
b - s1*rho
c - s2
### new:
f - function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c, eps0 = eps)
# new (04-25-06): eps0 instead of eps as local variable to f
{
nen0 - m2+c0*u
#catch a division by 0
nen - ifelse(abs(nen0)eps0, nen0, ifelse(nen00, nen0+eps0, nen0-eps0))
dnorm(u)/a0/nen * dnorm( t/a0/nen -(m1+b0*u)/a0)
}
-integrate(f, -Inf, -m2/c, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value+
integrate(f, -m2/c, Inf, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value
}
###
#
#cdf of X1X2, X1~N(m1,s1^2), X2~N(m2,s2^2), corr(X1,X2)=rho, evaluated at t
#
# eps is a very small number to catch errors in division by 0
###
#
pnnorm - function(t, m1, m2, s1, s2, rho, eps = .Machine$double.eps ^ 0.5){
a - s1*sqrt(1-rho^2)
b - s1*rho
c - s2
### new:
fp - function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c, eps0 = eps)
# new (04-25-06): eps0 instead of eps as local variable to fp
{nen0 - m2+c0*u ## for all u's used in integrate: never negative
#catch a division by 0
nen - ifelse(nen0eps0, nen0, nen0+eps0)
dnorm(u) * pnorm( t/a0/nen- (m1+b0*u)/a0)
}
### new:
fm - function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c, eps0 = eps)
# new (04-25-06): eps0 instead of eps as local variable to fm
{nen0 - m2+c0*u ## for all u's used in integrate: never positive
#catch a division by 0
nen - ifelse(nen0 (-eps0), nen0, nen0-eps0)
dnorm(u) * pnorm(-t/a0/nen+ (m1+b0*u)/a0)
}
integrate(fm, -Inf, -m2/c, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value+
integrate(fp, -m2/c, Inf, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value
}
##
For me this gives, e.g.:
pnnorm(0.5,m1=2,m2=3,s1=2,s2=1.4,rho=0.8)
[1] 0.1891655
dnnorm(0.5,m1=2,m2=3,s1=2,s2=1.4,rho=0.8)
[1] 0.07805282
Hth, Peter
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] distribution of the product of two correlated normal
Yu, Xuesong writes:
Does anyone know what the distribution for the product of two correlated
normal? Say I have X~N(a, \sigma1^2) and Y~N(b, \sigma2^2), and the
\rou(X,Y) is not equal to 0, I want to know the pdf or cdf of XY. Thanks
a lot in advance.
There is no closed-form expression (at least not to my knowledge) ---
but you could easily write some code for a numerical evaluation of the pdf /
cdf:
###
#code by P. Ruckdeschel, [EMAIL PROTECTED] 04-24-06
###
#
#pdf of X1X2, X1~N(m1,s1^2), X2~N(m2,s2^2), corr(X1,X2)=rho, evaluated at t
#
# eps is a very small number to catch errors in division by 0
###
#
dnnorm - function(t, m1, m2, s1, s2, rho, eps = .Machine$double.eps ^ 0.5){
a - s1*sqrt(1-rho^2)
b - s1*rho
c - s2
f - function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c, eps = eps)
{
nen0 - m2+c0*u
#catch a division by 0
nen - ifelse(abs(nen0)eps, nen0, ifelse(nen00, nen0+eps, nen0-eps))
dnorm(u)/a0/nen * dnorm( t/a0/nen -(m1+b0*u)/a0)
}
-integrate(f, -Inf, -m2/c, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value+
integrate(f, -m2/c, Inf, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value
}
###
#
#cdf of X1X2, X1~N(m1,s1^2), X2~N(m2,s2^2), corr(X1,X2)=rho, evaluated at t
#
# eps is a very small number to catch errors in division by 0
###
#
pnnorm - function(t, m1, m2, s1, s2, rho, eps = .Machine$double.eps ^ 0.5){
a - s1*sqrt(1-rho^2)
b - s1*rho
c - s2
fp - function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c, eps = eps)
{nen0 - m2+c0*u ## for all u's used in integrate: never negative
#catch a division by 0
nen - ifelse(nen0eps, nen0, nen0+eps)
dnorm(u) * pnorm( t/a0/nen- (m1+b0*u)/a0)
}
fm - function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c, eps = eps)
{
nen0 - m2+c0*u ## for all u's used in integrate: never positive
#catch a division by 0
nen - ifelse(nen0 -eps, nen0, nen0-eps)
dnorm(u) * pnorm(-t/a0/nen+ (m1+b0*u)/a0)
}
integrate(fm, -Inf, -m2/c, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value+
integrate(fp, -m2/c, Inf, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value
}
##
If you have to evalute dnnorm() or pnnorm() at a lot of values of t
for some given m1, m2, s1, s2, rho, then you should first evaluate
[p,d]nnorm() on a (smaller) number of gridpoints of values for t first
and then use something like approxfun() or splinefun() to give you a
much faster evaluable function.
Hth, Peter
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] distribution of the product of two correlated normal
Hi, Does anyone know what the distribution for the product of two correlated normal? Say I have X~N(a, \sigma1^2) and Y~N(b, \sigma2^2), and the \rou(X,Y) is not equal to 0, I want to know the pdf or cdf of XY. Thanks a lot in advance. yu [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
