«
»

Solve double integrals with R numerically

The goal of this post is to show how to evaluate in R terms like

\begin{displaymath }\int_a^b\int_c^d f(x,y)\; dx\; dy \end{displaymath}

with
\begin{displaymath}a,b,c,d\in[\infty,\infty]=:\overline{\ensuremath{\mathbb{R}}\...mathbb{R}}\xspace ^2\rightarrow\ensuremath{\mathbb{R}}\xspace .\end{displaymath}

The R function “integrate”:
We use the R function integrate and cite some arguments from its help:
integrate(f, lower, upper, …)


       f: an R function taking a numeric first argument and 
          returning a numeric vector of the same length.  
          Returning a non-finite element will generate an error.


lower, upper: the limits of integration.  Can be infinite.


     …: additional arguments to be passed to ‘f’.

There are further arguments, but see ?integrate for the complete and correct signature of this function.

Examples:

\begin{displaymath}\int_0^3x^2\;dx = 9\end{displaymath}

> integrate(function(x){x^2},0,3)
9 with absolute error < 1e-13

\begin{displaymath}\int_0^{2\pi}\sin(x)\;dx = 0\end{displaymath}

> integrate(sin,0,2*pi)
2.032977e-16 with absolute error < 4.4e-14


\begin{displaymath}\int_1^{\infty}\frac{1}{x^2}\;dx = 1\end{displaymath}

> integrate(function(x){1/x^2},1,Inf)
1 with absolute error < 1.1e-14

Double integral:
To evaluate
\begin{displaymath}\int_a^b\int_c^d f(x,y)\; dx\; dy\end{displaymath}

we define:
\begin{displaymath}g:\;y \mapsto \int_c^d f(x,y)\; dx\end{displaymath}

g <- function(y) {integrate(function(x) {f(x,y)},c,d)$value}

and then integrate over g:

\begin{displaymath}\int_a^b\int_c^d f(x,y)\; dx\; dy=\int_a^b g(y)\; dy\end{displaymath}
integrate(Vectorize(g),a,b)$value

Example:
$f(x,y)=x^2\cdot y,\;a=c=0,\;b=d=1$,
\begin{displaymath}\int_0^1\int_0^1 x^2\cdot y\; dx\; dy=...\end{displaymath}

R-Code:
> f <- function(x,y) { x^2*y }
> g <- (function(y) {integrate(function(x) 
{f(x,y)},0,1)$value})
> integrate(Vectorize(g),0,1)
0.1666667 with absolute error < 1.9e-15

Warning
Not all functions are suited: “Like all numerical integration routines, these evaluate the function on a finite set of points. If the function is approximately constant (in particular, zero) over nearly all its range it is possible that the result and error estimate may be seriously wrong” (from the integrate help – see ?integrate).

To demonstrate this, we integrate over the indicator function of ]1,3[:

plot

> f<-function(x){ifelse(x<1,0,ifelse(x<3,1,0))}
> integrate(f,-20,20)
2 with absolute error < 2.2e-15
> integrate(f,-30,30)
0 with absolute error < 0

On the other hand the last integral

\begin{displaymath}\int_0^\infty\left(\int_{-\infty}^\infty\binom{m}{j}[\Phi(b...\d t\right)f_\nu(s)\d s\end{displaymath}

we were interested in was no problem to evaluate… ;-)

Since LaTeX2HTML2Wordpress is a stressful conversion, I defer the theoretical background indefinitely.

Tags: , ,

Dieser Artikel wurde am Mittwoch, 06. Februar 2008 um 22:04 erstellt und ist in der Kategorie Computer abgelegt. Antworten zu diesem Artikel können durch den RSS 2.0-Feed verfolgt werden. Es besteht die Möglichkeit auf diesen Artikel zu antworten oder einen Trackback von der eigenen Seite zu senden.

3 Antworten zu “Solve double integrals with R numerically”

  1. Irmtraut sagt:
    9. Februar 2008 um 21:43
    Mozilla Firefox 1.5.0.2 Linux

    Was bedeutet das alles? Und für wen hast Du das geschrieben oder nur für Dich selbst? Wann kommt der nächste blogbeitrag für normale Sterbliche?

    Antworten
  2. felix sagt:
    10. Februar 2008 um 15:23
    Mozilla Firefox 2.0.0.9 Mac OS

    hey kornel, now you’re also posting in english? :)
    maybe this page is of interest to you: http://www-public.tu-bs.de:8080/~petras/teaching.html
    the lecture notes on “ Numerische Methoden in der Finanzmathematik” contain a chapter on numerical integration in several dimensions (with emphasis on really high dimension); this is not exactly what you are doing here, but maybe also interesting in this case. (he also gave this lecture in oldenburg, but it’s not listed on that page anymore…)
    the best approach seems to be to use (quasi) monte carlo simulations, at least in the case of high dimensions…
    ciao, felix

    Antworten
  3. Kornel sagt:
    13. Februar 2008 um 22:54
    Mozilla Firefox 2.0.0.12 Ubuntu Linux

    Hallo liebe Mutter.
    Für wen ich das geschrieben habe? Für Leute, die in R einfache Doppelintegrale lösen wollen. Es ist einfach, ja. Aber erstaunlicher Weise konnte das keiner bei mir auf der Arbeit.
    Du glaubst gar nicht wie viele Leute über Google hier auf die Seite kommen – einer hat schon nach “R Integrals” gesucht und ist hier gelandet.
    Und außerdem, ja, ist es eine Notiz für mich selber.

    Hi Felix.
    thank you for the link and comment. Actually I am interested in numerical integration in high dimensions. But only for some some special distributions and certainly not with R’s integrate function.
    I wanted to write more about that – but… when I have time and as a own post. Thanks for the patience with me. :-)
    Kornel.

    Antworten

Eine Antwort hinterlassen