I'm struggling on how I can simplify the quotient of two normal probability functions in R. Actually, I'm calculating a conditional skew-Normal density, them I have the division between this two function:
pnorm(alpha0+t(alpha2)%*%chol2inv(chol(omega2))%*%t(y2-xi2.1))/pnorm(tau2.1)
where alpha0+t(alpha2)%*%chol2inv(chol(omega2))%*%t(y2-xi2.1) and tau2.1 result in real numbers. For example, sometimes I have pnorm(-50)/pnorm(-40), e.g. an inconsistency 0/0. But these values are not zero, R is just approximating. I tried to use the erf function, but I got the same problem (0/0).
Any hint on how can I overcome this issue?
pnorm has a log parameter, which makes it return log(p). Change your equation to exp(log(p1) - log(p2)):
exp(pnorm(-50, log = TRUE) - pnorm(-40, log = TRUE))
#[1] 2.95577e-196