DeaR list,
I'm running an external program that computes some electromagnetic
response of a scattering body. The numerical scheme is based on a
discretization with a characteristic mesh size "y". The smaller y is,
the better the result (but obviously the computation will take longer).
A convergence study showed the error between the computed values and
the exact solution of the problem to be a quadratic in y, with
standard error increasing as y^3. I wrote the interface to the
program in R, as it is much more user friendly and allows for post-
processing analysis. Currently, it only runs with user-defined
discretization parameter. I would like to implement an adaptive
scheme [1] and provide the following improvements,
1) obtain an estimate of the error by fitting the result against a
series of mesh sizes with the quadratic model, and extrapolate at y =
0. (quite straight forward)
2) adapt "dynamically" the set of mesh sizes to fulfill a final
accuracy condition, between a starting value (a rule-of thumb
estimate is given by the problem values). The lower limit of y should
also be constrained by the resources (again, an empirical rule
dictates the computation time and memory usage).
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@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.
