Think about what you are actually getting from pnorm().
For each scalar "s" you want to get
pnorm(s,mu[2],sigma[2]) pnorm(s,mu[3],sigma[3])
pnorm(s,mu[4],sigma[4])
(and then take products of things).
But when "s" is a vector you get
pnorm(s[1],mu[2],sigma[2]) pnorm(s[2],mu[3],sigma[3])
pnorm(s[3],mu[4],sigma[4])
pnorm(s[4],mu[2],sigma[2]) pnorm(s[5],mu[3],sigma[3])
pnorm(s[6],mu[4],sigma[4]) ...
and the product will be nothing like what you want.
Vectors get re-cycled in R.
Note that you get a single scalar quantity from the "prod(1 -
pnorm(...))" component of
test2(). That single scalar multiplies each entry of the vector
dnorm(s,mu[1],sigma[1]).
Whereas you want dnorm(s[i],mu[1],sigma[1]) to be multiplied by a
product involving
only pnorm() terms evaluated at s[i]. That is what Vectorize() arranges
for you.
If you are going to use R you should learn its syntax and know about
things like the
re-cycling of vectors.
cheers,
Rolf Turner
On 11/22/13 04:46, dan wang wrote:
Hi all,
I tried below two methods to calculate the integrate of a same function.
Two different results are given.
The first one should be the right answer.
Can any one help to explain why?
Another issue is, the first one is not convenient as I have to update the
mu and sigma outside the function. How can I modify the second one so I
could issue the parameters in the integrate function.
Thanks in advance,
Dan
test <- function(s){
prod(1-pnorm(s,mu[-1],sigma[-1]))*dnorm(s,mu[1],sigma[1])
}
testV <- Vectorize(test)
mu=c(0,0,0,0)
sigma=c(1,1,1,1)
integrate(testV,lower=-Inf,upper=Inf)$value
test2 <- function(s, mu, sigma){
prod(1-pnorm(s,mu[-1],sigma[-1]))*dnorm(s,mu[1],sigma[1])
}
integrate(test2,mu=c(0,0,0,0),sigma=c(1,1,1,1),lower=-Inf,upper=Inf)$value
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
