[R] Numerical integration

[Dennis Murphy]

Hi:

You could write the function this way:

f <- function(x, xstar, k) dnorm(x) * k * x * (x >= xstar)

where the term in parentheses is a logical. For any x < xstar, f will
be zero by definition. Substitute your density in for dnorm().

To integrate over a grid of (xstar, k) values, you could try this:

# generate a grid of (xstar, k) pairs
pars <- expand.grid(xstar = seq(0, 2, by = 0.2), k = seq(0.5, 2, by = 0.5))

# We modify the function to do the integration and return its value.
# Note that since f = 0 for all x <= xstar, integrating from -Inf to Inf
# yields the same answer as integrating from xstar to Inf.
fun <- function(xstar, k) {
    f <- function(x, xstar, k) dnorm(x) * k * x * (x >= xstar)
    integrate(f, -Inf, Inf, xstar = xstar, k = k)$value
   }

# Method 1: (outputs a data frame)
library(plyr)
out <- mdply(pars, fun)

# Method 2a: (outputs a matrix)
out <- cbind(pars, mapply(fun, xstar = pars[['xstar']], k = pars[['k']]))

# Method 2b: (outputs a data frame)
out <- data.frame(pars, value = mapply(fun, xstar = pars[['xstar']], k
= pars[['k']]))

HTH,
Dennis


On Wed, Jun 29, 2011 at 4:52 PM, nany23 <anna.botto at gmail.com> wrote:
> Hello!
>
> I know that probably my question is rather simple  but I' m a very beginner
> R-user.
>
> I have to numerically integrate  the product of two function A(x) and B(x).
> The integretion limits are [X*; +inf]
>
> Function A(x) is a pdf function while B(x)=e*x is a linear function whose
> value is equal to 0 when the  x < X*
>
> Moreover I have to iterate this process for different value of X* and for
> different pdf of the same type.
>
> I know the comand INTEGRATE but  I can' t make it work.
>
> Which is the best function to do this and how does it work?
>
> Thank you very much in advanced!!
>
> --
> View this message in context: http://r.789695.n4.nabble.com/Numerical-integration-tp3634365p3634365.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at 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.
>