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)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 181 bytes
Desc: Digital signature
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20140801/18ace718/attachment.bin>