[R] Numerical integration with R
David Winsemius
dwinsemius at comcast.net
Sat Jun 14 22:18:22 CEST 2014
On Jun 14, 2014, at 7:40 AM, Christofer Bogaso wrote:
> Hi again,
>
> I was tying to solve following 2-fold integration with package cubature.
> However spending approximately 2 hours it failed to generate any number. I
> am using latest R with win-7 machine having 4gb ram.
>
>> library(cubature)
>> f <- function(x) {
> + z1 <- x[1]
> + z2 <- x[2]
> +
> + Rho = 1
> +
> + L <- log(1 - pnorm(z1)) * log(1 - pnorm(Rho * z1 + sqrt(1 - Rho^2) * z2))
> * dnorm(z1) * dnorm(z2)
> +
> + return(L)
> + }; adaptIntegrate(f, lowerLimit = c(-Inf, -Inf), upperLimit = c(Inf, Inf))
> Error in adaptIntegrate(f, lowerLimit = c(-Inf, -Inf), upperLimit = c(Inf,
> :
> 'Realloc' could not re-allocate memory (1006632880 bytes)
> In addition: Warning messages:
> 1: In adaptIntegrate(f, lowerLimit = c(-Inf, -Inf), upperLimit = c(Inf, :
> Reached total allocation of 1535Mb: see help(memory.size)
> 2: In adaptIntegrate(f, lowerLimit = c(-Inf, -Inf), upperLimit = c(Inf, :
> Reached total allocation of 1535Mb: see help(memory.size)
>
>
> Would really appreciate if someone helps me to sort out this problem.
You are posting a problem with no data. You are also not explaining what you are trying to do. In particular, what was your intent when you wrote:
upperLimit = c(Inf, Inf) # ????
When I run this with finite limits that are increasing in range and several values for a length 2 "x"-vector
x <- c(0.1,0.1)
x <-c(.1, 0 )
x <- c(10,10)
x <- c(10, -10)
I consistently get:
#------------
> adaptIntegrate(f, lowerLimit = c(-1000, 1000), upperLimit = c(-1000, 1000))
$integral
[1] 0
$error
[1] 0
$functionEvaluations
[1] 17
$returnCode
[1] 0
#---------------
Which makes me think there is something mathematically degenerate about your "L" expression.
>
> Thanks and regards,
>
> [[alternative HTML version deleted]]
Arrrgh. Pray tell, exactly how many years have you been posting to R-help? I've checked... and it's really not that difficult to get gmail to send plain text.
--
David Winsemius
Alameda, CA, USA
More information about the R-help
mailing list