4

It's my first time on stackoverflow, so please be patient with me! I'm using R-package "pracma" to calculate a triple integral. This is my code:

mu1=0
mu2=0
mu3=0
mu4=0
sigma1=1
sigma2=1
sigma3=1
sigma4=1
f3=function(x,y,z){dnorm(z,mean = mu2,sd=sigma2)*dnorm(y,mean = 
mu3,sd=sigma3)*(1-pnorm(z,mean= mu1,sd=sigma1))*dnorm(x,mean = 
mu4,sd=sigma4)}
ymin=function(x){x}
zmin=function(x,y){y}
integral3(f3,xmin = -100,xmax = 100,ymin,ymax=100,zmin,zmax = 100)

But, actually, I obtain the error below:

 Error in if (adjerr[1] > localtol) { :

  missing value where TRUE/FALSE needed

Can someone give me a cue to solve this problem? Thank in advance

Elena
  • 41
  • 1

1 Answers1

0

This happens because of your max and min values you set. In your f3() function you operate with dnorm() which reaches 0 very quickly once you are outside [-8,8] interval:

dnorm(0:20)
##[1] 3.989423e-01 2.419707e-01 5.399097e-02 4.431848e-03 1.338302e-04 1.486720e-06 6.075883e-09 9.134720e-12 5.052271e-15 1.027977e-18
##[11] 7.694599e-23 2.118819e-27 2.146384e-32 7.998828e-38 1.096607e-43 5.530710e-50 1.026163e-56 7.004182e-64 1.758750e-71 1.624636e-79
##[21] 5.520948e-88

If you use some reasonable values for your min and max parameters, the function works as expected:

library(pracma)
mu1=0
mu2=0
mu3=0
mu4=0
sigma1=1
sigma2=1
sigma3=1
sigma4=1
f3=function(x,y,z){dnorm(z,mean = mu2,sd=sigma2)*dnorm(y,mean = 
                                                         mu3,sd=sigma3)*(1-pnorm(z,mean= mu1,sd=sigma1))*dnorm(x,mean = 
                                                                                                                 mu4,sd=sigma4)}
xmin <- -10
xmax <-  10
ymin=function(x){x}
ymax <- 10
zmin=function(x,y){y}
zmax <- 10
integral3(f3, xmin, xmax ,ymin, xmax, zmin, zmax)
# [1] 0.04166667
Katia
  • 3,784
  • 1
  • 14
  • 27
  • Thanks a lot! I set those max and min values because my intention was to integrate between minus infinite and infinite. – Elena May 18 '18 at 16:22
  • @Elena Some internal functions that are used in this package specifically .tensor() and .save2list() calculate values that will not evaluate and will cause this error. You can look at the source code here: https://github.com/cran/pracma/blob/master/R/integral2.R (integral2() is called from within integral3()). – Katia May 18 '18 at 16:37