Michael Lindner
2007-Aug-23 06:07 UTC
[R] nls() and numerical integration (e.g. integrate()) working together?
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
Frede Aakmann Tøgersen
2007-Aug-23 07:05 UTC
[R] nls() and numerical integration (e.g. integrate()) workingtogether?
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.
>