CODEWITHSUNDEEP
rmathstatistics
David
David

Reputation: 213

Problems with Gaussian Quadrature in R

I'm using the the gaussquad package to evaluate some integrals numerically.

I thought the ghermite.h.quadrature command worked by evaluating a function f(x) at points x1, ..., xn and then constructing the sum w1*f(x1) + ... + wn*f(xn), where x1, ..., xn and w1, ..., wn are nodes and weights supplied by the user.

Thus I thought the commands

ghermite.h.quadrature(f,rule)
sum(sapply(rule$x,f)*rule$w)

would yield the same output for any function f, where ''rule'' is a dataframe which stores the nodes in a column labelled ''x'' and the weights in a column labelled "w". For many functions the output is indeed the same, but for some functions I get very different results. Can someone please help me understand this discrepancy?

Thanks!

Code:

n.quad = 50 
rule = hermite.h.quadrature.rules(n.quad)[[n.quad]]

f <- function(z){ 

 f1 <- function(x,y) pnorm(x+y)

 f2 <- function(y) ghermite.h.quadrature(f1,rule,y = y)

 g <- function(x,y) x/(1+y) / f2(y)*pnorm(x+y)

 h <- function(y) ghermite.h.quadrature(g,rule,y=y)

 h(z)

}

ghermite.h.quadrature(f,rule)
sum(sapply(rule$x,f)*rule$w)

Upvotes: 1

Views: 1137

Answers (1)

Severin Pappadeux
Severin Pappadeux

Reputation: 20080

Ok, that problem got me interested.

I've looked into gaussquad sources, and clearly author is not running sapply internally, because all integrands/function shall return vector on vector argument.

It is clearly stated in documentation:

functn an R function which should take a numeric argument x and possibly some parameters. The function returns a numerical vector value for the given argument x

In case where you're using some internal functions, they're written that way, so everything works.

You have to rewrite your function to work with vector argument and return back a vector

UPDATE

Vectorize() works for me to rectify the problem, as well as simple wrapper with sapply

vf <- function(z) {
    sapply(z, f)
}

After either of those changes, results are identical: 0.2029512

Upvotes: 1

Related Questions