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