R 敏感性分析

雖然在R中已經有sensitivity這個套件可以用來分析了，

但還是動手做出自己想要的比較有成就感。

首先定一個要測試的函數

fit <- function(x) (x[2] - x[1]^2)^2 + (1 - x[1])^2

設定有哪些自變數及變動的範圍

temp <- list(a = seq(-2, 4, 0.5), b = seq(-1, 4, length.out = 11))
temp

## $a
##  [1] -2.0 -1.5 -1.0 -0.5  0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0
## 
## $b
##  [1] -1.0 -0.5  0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0
## 

先產生一個List（自變數的個數×自變數的個數×變動幅度的長度）

temp2 <- lapply(temp, function(x) matrix(rep(sapply(temp, median), 
    length(x)), length(x), byrow = T))
temp2

## $a
##       [,1] [,2]
##  [1,]    1  1.5
##  [2,]    1  1.5
##  [3,]    1  1.5
##  [4,]    1  1.5
##  [5,]    1  1.5
##  [6,]    1  1.5
##  [7,]    1  1.5
##  [8,]    1  1.5
##  [9,]    1  1.5
## [10,]    1  1.5
## [11,]    1  1.5
## [12,]    1  1.5
## [13,]    1  1.5
## 
## $b
##       [,1] [,2]
##  [1,]    1  1.5
##  [2,]    1  1.5
##  [3,]    1  1.5
##  [4,]    1  1.5
##  [5,]    1  1.5
##  [6,]    1  1.5
##  [7,]    1  1.5
##  [8,]    1  1.5
##  [9,]    1  1.5
## [10,]    1  1.5
## [11,]    1  1.5
## 

把要測量的變動幅度帶進去List

for (i in 1:length(temp)) {
    temp2[[i]][, i] <- temp[[i]]
}
temp2

## $a
##       [,1] [,2]
##  [1,] -2.0  1.5
##  [2,] -1.5  1.5
##  [3,] -1.0  1.5
##  [4,] -0.5  1.5
##  [5,]  0.0  1.5
##  [6,]  0.5  1.5
##  [7,]  1.0  1.5
##  [8,]  1.5  1.5
##  [9,]  2.0  1.5
## [10,]  2.5  1.5
## [11,]  3.0  1.5
## [12,]  3.5  1.5
## [13,]  4.0  1.5
## 
## $b
##       [,1] [,2]
##  [1,]    1 -1.0
##  [2,]    1 -0.5
##  [3,]    1  0.0
##  [4,]    1  0.5
##  [5,]    1  1.0
##  [6,]    1  1.5
##  [7,]    1  2.0
##  [8,]    1  2.5
##  [9,]    1  3.0
## [10,]    1  3.5
## [11,]    1  4.0
## 

計算出每個自變數每次變動的值，作為Y座標

ytemp <- lapply(temp2, function(x) apply(x, 1, fit))
ytemp

## $a
##  [1]  15.2500   6.8125   4.2500   3.8125   3.2500   1.8125   0.2500
##  [8]   0.8125   7.2500  24.8125  60.2500 121.8125 219.2500
## 
## $b
##  [1] 4.00 2.25 1.00 0.25 0.00 0.25 1.00 2.25 4.00 6.25 9.00
## 

計算變動的百分比，作為X座標

xtemp <- lapply(temp, function(x) (x - median(x))/median(x))
xtemp

## $a
##  [1] -3.0 -2.5 -2.0 -1.5 -1.0 -0.5  0.0  0.5  1.0  1.5  2.0  2.5  3.0
## 
## $b
##  [1] -1.6667 -1.3333 -1.0000 -0.6667 -0.3333  0.0000  0.3333  0.6667
##  [9]  1.0000  1.3333  1.6667
## 

先畫出第一條自變數的線

plot(xtemp[[1]], ytemp[[1]], type = "o", xlab = "", ylab = "")

再畫出其他自變數的線

mapply(function(x) {
    points(xtemp[[x]], ytemp[[x]], pch = x)
    points(xtemp[[x]], ytemp[[x]], type = "l")
}, 2:length(temp))

再加上個說明圖

legend(0, max(unlist(ytemp)), names(temp), pch = c(1:length(temp)))

這樣簡易的敏感性分析就完成啦！