[Sensitivity-commits] r21 - pkg/R

[noreply at r-forge.r-project.org]

Author: pujol
Date: 2010-02-10 18:55:27 +0100 (Wed, 10 Feb 2010)
New Revision: 21

Modified:
   pkg/R/morris.R
   pkg/R/morris_oat.R
Log:
Some bug corrections (Morris OAT) :
* design randomization
* scaling : scale=FALSE by def, no more y scaling, and X scaling by homothety

Modified: pkg/R/morris.R
===================================================================
--- pkg/R/morris.R	2010-02-08 20:02:25 UTC (rev 20)
+++ pkg/R/morris.R	2010-02-10 17:55:27 UTC (rev 21)
@@ -1,8 +1,6 @@
-# Morris's screening method (Morris 1992, Campolongo 2007)
-# Provides also simplex-based screening designs (Pujol 2007)
+# Morris's screening method (M. D. Morris 1992, F. Campolongo 2007)
+# Provides also simplex-based screening designs (G. Pujol 2007)
 #
-# Gilles Pujol 2006-2008
-#
 # Sub-files:
 # * morris_oat.R
 # * morris_simplex.R
@@ -10,13 +8,13 @@
 
 
 ind.rep <- function(i, p) {
-# indices of the ith trajectory in the DoE
+# indices of the points of the ith trajectory in the DoE
   (1 : (p + 1)) + (i - 1) * (p + 1)
 }
 
 
 morris <- function(model = NULL, factors, r, design, binf = 0, bsup = 1,
-                   scale = TRUE, ...) {
+                   scale = FALSE, ...) {
   
   # argument checking: factor number and names
   if (is.character(factors)) {
@@ -126,8 +124,11 @@
   y <- x$y
 
   if (x$scale) {
-    X <- scale(X)
-    y <- as.numeric(scale(y))
+    #X <- scale(X)
+	#y <- as.numeric(scale(y))
+	Binf <- matrix(x$binf, nrow = nrow(X), ncol = length(x$binf), byrow = TRUE)
+	Bsup <- matrix(x$bsup, nrow = nrow(X), ncol = length(x$bsup), byrow = TRUE)
+	X <- (X - Binf) / (Bsup - Binf) 
   }
 
   if (x$design$type == "oat") {

Modified: pkg/R/morris_oat.R
===================================================================
--- pkg/R/morris_oat.R	2010-02-08 20:02:25 UTC (rev 20)
+++ pkg/R/morris_oat.R	2010-02-10 17:55:27 UTC (rev 21)
@@ -1,18 +1,25 @@
-# Morris' OAT sub-routines. (See main file morris.R)
-#
-# Gilles Pujol 2006
+# Morris's OAT method sub-routines. (See main file morris.R)
 
-
-random.oat <- function(p, nl, jump) {
-  delta <- jump / (nl - 1) 
-  B <- matrix(0, nrow = p + 1, ncol = p)                 # orientation matrix
+random.oat <- function(p, nl, design.step) {
+  # orientation matrix B
+  B <- matrix(-1, nrow = p + 1, ncol = p)                 
   B[lower.tri(B)] <- 1
-  x.base <- matrix(nrow = p + 1, ncol = p)               # base point
-  D <- diag(sample(c(-1, 1), size = p, replace = TRUE))  # directions 
+  # directions matrix D
+  D <- diag(sample(c(-1, 1), size = p, replace = TRUE)) 
+  # permutation matrix P
+  perm <- sample(p)
+  P <- matrix(0, nrow = p, ncol = p)
+  for (i in 1 : p) {
+    P[i, perm[i]] <- 1
+  }
+  # starting point
+  x.base <- matrix(nrow = p + 1, ncol = p)               
   for (i in 1 : p) {                                    
-    x.base[,i] <- ((sample(nl[i] - jump[i], size = 1) - 1) / (nl[i] - 1))
+    x.base[,i] <- ((sample(nl[i] - design.step[i], size = 1) - 1) / (nl[i] - 1))
   }
-  0.5 * ((2 * B - 1) %*% D + 1) %*% diag(delta) + x.base
+  # grid step
+  delta <- design.step / (nl - 1)
+  return(0.5 * (B %*% P %*% D + 1) %*% diag(delta) + x.base)
 }
 
 
@@ -26,7 +33,9 @@
     j <- ind.rep(i, p)
     j1 <- j[1 : p]
     j2 <- j[2 : (p + 1)]
-    ee[i,] <- (y[j2] - y[j1]) / rowSums(X[j2,] - X[j1,])
+	A <- X[j2,] - X[j1,]
+	b <- y[j2] - y[j1]
+	ee[i,] <- solve(A, b)
   }
   return(ee)
 }