From: baptiste Auguié <ba208_at_exeter.ac.uk>

Date: Sun, 04 May 2008 16:39:43 +0100

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. Received on Sun 04 May 2008 - 16:38:59 GMT

Date: Sun, 04 May 2008 16:39:43 +0100

I'm looking for advice on this second point (both on the technical aspect, and whether this is sound statistically):

- I can foresee that I should always start with a few y values before I can do any extrapolation, but how many of them? 3, 10? How could I know?
- once I have enough points (say, 10) to use the fitting procedure and get an estimate of the error, how should I decide the best location of the next y if the error is too important?
- in a practical implementation, I would use a while loop and append the successive values to a data.frame(y, value). However, this procedure will be run for different parameters (wavelengths, actually), so the set and number of y values may vary between one run and another. I think I'd be better off using a list with each new run having its own data.frame. Does this make sense?

Below are a few lines of code to illustrate the problem,

*> program.result <- function(x, p){ # made up function that mimicks
**> the results of the real program
**>
**> y <- p[3]*x^2 + p[2]*x + p[1]
**> y * (1 + rnorm(1, mean=0, sd = 0.1 * y^3))
**> }
**>
**>
**> p0 <- c(0.1, 0.1, 2) # set of parameters
**>
**> ## user defined limits of the y parameter (log scale)
**> limits <- c(0.1, 0.8)
**> limits.log <- (10^limits)
**> y.log <- seq(limits.log[1], limits.log[2], l=10)
**>
**> y <- log10(y.log)
**>
**> result <- sapply(y, function(x) program.result(x, p0)) # results of
**> the program
**>
**> #### fitting and extrapolation procedure ####
**> library(gplots) # plot with CI
**> plotCI(y, result, y^3, xlim=c(0, 1), ylim=c(0, 2)) # the data
**> with y^3 errors
**>
**> my.data <- data.frame(y = y, value = result)
**>
**> fm <- lm(value ~ poly(y, degree=2, raw=TRUE), data = my.data ,
**> weights = 1/y^3)
**> lines(y, predict(fm, data.frame(y=y)), col = 2)
**>
**> extrap <- summary(fm)$coefficients[1,] # intercept and error on it
**> plotCI(0,extrap[1], 2 * extrap[2], col = 2, add=T)
**>
**>
**> ### my naive take on adaptive runs... ##
**>
**> objective <- 1e-3 # stop when the standard error of the
**> extrapolated value is smaller than this
**>
**> err <- extrap[2]
**> my.color <- 3
**> while (err > objective){
**>
**> new.value <- min(y)/2 # i don't know how to choose this optimally
**> y <- c(new.value, y)
**> new.result <- program.result(new.value, p0)
**> result <- c(new.result, result)
**> points(new.value, new.result, col= my.color)
**> my.data <- data.frame(y = y, value = result)
**> fm <- lm(value ~ poly(y, degree=2, raw=TRUE), data = my.data ,
**> weights = 1/y^3)
**> lines(y, predict(fm, data.frame(y=y)), col = my.color)
**>
**> extrap <- summary(fm)$coefficients[1,] # intercept and error on it
**> err <- extrap[2]
**> print(err)
**> plotCI(0,extrap[1], 2 * err, col = 2, add=T)
**> my.color <- my.color + 1
**>
**> }
**> err
*

Many thanks in advance for your comments,

baptiste

[1]: Yurkin et al., Convergence of the discrete dipole approximation. II. An extrapolation technique to increase the accuracy. J. Opt. Soc. Am. A / Vol. 23, No. 10 / October 2006

Baptiste Auguié

Physics Department

University of Exeter

Stocker Road,

Exeter, Devon,

EX4 4QL, UK

Phone: +44 1392 264187

http://newton.ex.ac.uk/research/emag

http://projects.ex.ac.uk/atto

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. Received on Sun 04 May 2008 - 16:38:59 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Tue 06 May 2008 - 10:30:33 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*