Hi, I have been using the R/fast99 library to work through the Extended FAST paper by Saltelli et al 1999 [1]. In trying to repeat some examples with R, I have encountered the following 'bug' in the implementation.
Synopsis: tell.fast99 gives NA first order sensitivity indices when 'M *
omega[1] > (n / 2) - 1'.
This just so happens for n=65,M=4,omega[1]=8, which is described in the paper.
Offending tell.fast99 lines
===========================
for (i in 1 : p) {
l <- seq((i - 1) * n + 1, i * n)
f <- fft(x$y[l], inverse = FALSE)
Sp <- ( Mod(f[2 : (n / 2)]) / n )^2
^^^^^^^^^^^^^| SP length is (n / 2) - 1
V[i] <- 2 * sum(Sp)
D1[i] <- 2 * sum(Sp[(1 : x$M) * x$omega[1]])
^^^^^^^^^^^^^^^^^^^^^^^| now try and read SP[M *
omega[1]]
Dt[i] <- 2 * sum(Sp[1 : (x$omega[1] / 2)])
}
Worked example (copied from fast99 manpage, with n=65)
=======================================================
> x=fast99(model=ishigami.fun, n=65, factors=3,q="qunif",
> q.arg=list(min=-pi,max=pi))
> x$M
[1] 4
> x$omega
[1] 8 1 1
> x
Call:
fast99(model = ishigami.fun, factors = 3, n = 65, q = "qunif", q.arg = list(min
= -pi, max = pi))
Model runs: 195
Estimations of the indices:
first order total order
X1 NA 0.026476165
X2 NA 0.971861619
X3 NA 0.001241782
SimLab (the program that Saltelli uses in his books) gives consistent results
for n=65.
Why is the fft result truncated to '2 : (n / 2)'?
If there is a good reason, can it go in the manpage!
[1] A. Saltelli, S. Tarantola and K. Chan, 1999, A quantitative, model
independent method
for global sensitivity analysis of model output, Technometrics, 41, 39-56.
-- Peter
(A907 E02F A6E5 0CD2 34CD 20D2 6760 79C5 AC40 DD6B)
signature.asc
Description: Digital signature
______________________________________________ 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.
