[R] R sensitivity/fast99 n=65

Peter Ross pross at xvid.org
Fri Aug 1 14:12:22 CEST 2014


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>


More information about the R-help mailing list