[R] nls() and numerical integration (e.g. integrate()) workingtogether?

Thu Aug 23 09:05:50 CEST 2007

Michael, it should be possible to get nls and integrate to work together. Buuuut, there are several problems you need to consider first. The most important is the definition of e.g. your f-function and how integrate() work. The easiest way to show you the problem is like in the following.

The integrate() has an ... argument through which you can pass additional arguments to your integrand. Let's define

f <- function(z,x1,b1) {z + b1*x1}

where z is the variable of integration and x1 and b1 is thought of as parameters. We pass the parameters x1 and b1 to f through the ... argument of integrate() later on (see ?integrate).

Now since integrate() is an adaptive integration rutine you'll find, that integrate() based on some rule decides to partition the interval of integration (here (0,1)) into a number of subintervals. integrate() then evaluates f() in the points defining the partition and uses some quadrature rule to estimate the integral and its associated error. If the error is above the tolerance then integrate() makes further subdivision. The most important consequence of this is that you do not know the length of the vector z that is passed to f() by integrate(). Most likely this length is different from the length of x1. You can see this by trying out the following code:

f <- function(z,x1,b1) {
  print(z)
  print(b1)
  print(x1)
  z + b1*x1
}

# f is a vectorized funtion: try f(c(2,3),x1,b1=.5), so why use Vectorize()?

Now integrate f from 0 to 1.

integrate(f,0,1,x1=x1,b1=.5)

Note that we have the parameters x1 and b1 to f through the ... argument. 

You'll see some error message relating to the problem with different lengths of vectors.

You need to resolve this problem before continuing.

Med venlig hilsen / Regards

Frede Aakmann Tøgersen
Forsker / Scientist

   AARHUS UNIVERSITET / UNIVERSITY OF AARHUS 
Det Jordbrugsvidenskabelige Fakultet / Faculty of Agricultural Sciences 
Inst. for Genetik og Bioteknologi / Dept. of Genetics and Biotechnology 
Blichers Allé 20, P.O. BOX 50 
DK-8830 Tjele 

Tel: +45 8999 1900 
Direct: +45 8999 1878 
Mobile: +45  
E-mail: FredeA.Togersen at agrsci.dk 
Web: www.agrsci.dk 

--------------------------------------------------------------------------------
Tilmeld dig DJF's nyhedsbrev / Subscribe Faculty of Agricultural Sciences Newsletter. 

Denne email kan indeholde fortrolig information. Enhver brug eller offentliggørelse af denne email uden skriftlig tilladelse fra DJF er ikke tilladt. Hvis De ikke er den tiltænkte adressat, bedes De venligst straks underrette DJF samt slette emailen.

This email may contain information that is confidential. Any use or publication of this email without written permission from Faculty of Agricultural Sciences is not allowed. If you are not the intended recipient, please notify Faculty of Agricultural Sciences immediately and delete this email.

> -----Oprindelig meddelelse-----
> Fra: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] På vegne af Michael Lindner
> Sendt: 23. august 2007 08:07
> Til: r-help at stat.math.ethz.ch
> Emne: [R] nls() and numerical integration (e.g. integrate()) 
> workingtogether?
> 
> Dear List-Members,
> 
> since 3 weeks I have been heavily working on reproducing the 
> results of an economic paper. The method there uses the 
> numerical solution of an integral within nonlinear least 
> squares. Within the integrand there is also some parameter to 
> estimate. Is that in the end possible to implement in R 
> [Originally it was done in GAUSS]? I'm nearly into giving up.
> 
> I constucted an example to showing the problems I face.
> 
> I have three questions - related to three errors shown below: 
> 1) How to make clear that in the integrand z is the 
> integration variable and
> b1 is a parameter while x1 is a data variable
> 2) and 3) How to set up a correct estimation of the integral?
> 
> library(stats)
> y <- c(2,15,24,21,5,6,)
> x1 <- c(2.21,5,3,5,2,1)
> x2 <- c(4.51,6,2,11,0.4,3)
> f <- function(z) {z + b1*x1}
> vf <- Vectorize(f)
> g <- function(z) {z + x1}
> vg <- Vectorize(f)
> 
> Error 1: 
> > nls(y ~ integrate(vf,0,1)+b2*x2,start=list(b1=0.5,b2=2))
> Error in function (z)  : object "b1" not found
> 
> Error 2:
> > nls(y ~ integrate(vg,0,1)+b2*x2,start=list(b1=0.5,b2=2))
> Error in integrate(vg, 0, 1) : REAL() can only be applied to 
> a 'numeric', not a 'list'
> 
> Error 3:
> > nls(y ~ integrate(g,0,1)+b2*x2,start=list(b1=0.5,b2=2))
> Error in integrate(g, 0, 1) + b2 * x2 : non-numeric argument 
> to binary operator In addition: Warning messages:
> 1: longer object length
>         is not a multiple of shorter object length in: z + x1
> 
> With a lot of thanks in advance,
> 
> Michael
> 
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>