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(nen0>0, 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(nen0>eps, 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