汎用最適化関数 nlm() の使用法

統計学では、最尤推定や最小自乗法等、ある目的関数の最適化を行なう必要が しばしば起こる。関数 nlm (Non-Linear Minimization) はニュートン法に基づく非線形関数の最小化関数で、r-help 記事によればしばしば optim() 関数の諸手法よりも良い結果を与えるらしい。


nlm() の使用法

命令書式

 nlm(f, p, hessian = FALSE, typsize=rep(1, length(p)), fscale=1,
        print.level = 0, ndigit=12, gradtol = 1e-6,
        stepmax = max(1000 * sqrt(sum((p/typsize)^2)), 1000),
        steptol = 1e-6, iterlim = 100, check.analyticals = TRUE, ...)

引数:

詳細

もし gradient や hessian が与えられても好ましくないモードや長さを与えるようなら、check.analyticals = TRUE (既定動作) が指定されていれば警告とともに無視される。 hessian は、 gradient が与えられていなかったり、チェックをパスしない限り、無視される。

原典で与えられている三種類の手法のうち、直線探索である手法 1 だけが使われる。

返り値

次の成分を含むリスト(str 関数を使うと簡潔に一覧できる)

例1 demo(nlm) より

Helical Valley Function. Page 362 Dennis + Schnabel

nlm1 <- function () {
  theta <- function(x1,x2) atan(x2,x1)/(2*pi)
  f <- function(x) {                          # 目的関数定義
    f1 <- 10*(x[3] - 10*theta(x[1],x[2]))
    f2 <- 10*(sqrt(x[1]^2+x[2]^2)-1)
    f3 <- x[3]
    return(f1^2+f2^2+f3^2)
  }
  # x3=0 における曲面を見る
  x <- seq(-1, 2, len=50)
  y <- seq(-1, 1, len=50)
  # z は x,y のすべての組み合わせに対する関数値を順に並べたベクトル
  z <- apply(as.matrix(expand.grid(x, y)), 1, function(x) f(c(x, 0))) 
  contour(x, y, matrix(log10(z), 50, 50))     # (対数値)曲面の等高線図
  nlm.f <- nlm(f, c(-1,0,0), hessian=TRUE, print=0)
  points(rbind(nlm.f$estim[1:2]), col = "red", pch = 20) # 等高線図に最適パラメータをプロット
  nlm.f   
}
> nlm.f <- nlm1()
> nlm.f            # 非線形最小化の結果 
$minimum           # 最小値
[1] 1.238237e-14 
$estimate          # その時のパラメータ値   
[1]  1.000000e-00  3.071089e-09 -6.063445e-09
$gradient          # 収束点での偏微分ベクトル gradient (ゼロに近いはず)
[1] -3.755775e-07  3.485886e-06 -2.202374e-06
$hessian           # 収束点での二階導関数行列 hessian (正定値符合行列のはず)
              [,1]          [,2]          [,3]
[1,]  2.000000e+02   -0.04065842  9.774602e-07
[2,] -4.065842e-02  506.60589960 -3.183099e+02
[3,]  9.774602e-07 -318.30988572  2.020000e+02
$code              # 収束判定コード (おそらく成功)
[1] 2
$iterations        # 終了までの反復回数 
[1] 27

> str(nlm.f)       # 結果の簡潔な要約
List of 6
 $ minimum   : num 1.24e-14
 $ estimate  : num [1:3]  1.00e-00  3.07e-09 -6.06e-09
 $ gradient  : num [1:3] -3.76e-07  3.49e-06 -2.20e-06
 $ hessian   : num [1:3, 1:3]  2.00e+02 -4.07e-02  9.77e-07 -4.07e-02 5.07e+02 ...
 $ code      : int 2
 $ iterations: int 27
> eigen(nlm.f$hessian) 
$values                 # hessian の固有値はすべて正(正定値符合行列!)
[1] 707.173143 200.000000   1.432756
$vectors
              [,1]          [,2]         [,3]
[1,]  6.782635e-05  1.000000e-00 0.0001091525
[2,] -8.460531e-01 -8.043624e-07 0.5330985699
[3,]  5.330986e-01 -1.285070e-04 0.8460531324

#ref(): The style ref(filename,pagename) is ambiguous and become obsolete. Please try ref(pagename/filename)

例2 demo(nlm) より

### the Rosenbrock banana valley function

nlm2 <- function () {
  fR <- function(x) {  # 目的関数定義
    x1 <- x[1]; x2 <- x[2]
    100*(x2 - x1*x1)^2 + (1-x1)^2
  }
 fx <- function(x) {   # fR のベクトル化版 (z をスマートに定義するため)
   x1 <- x[,1]; x2 <- x[,2]
   100*(x2 - x1*x1)^2 + (1-x1)^2
 }
 x <- seq(-2, 2, len=100)
 y <- seq(-.5, 1.5, len=100)
 z <- fx(expand.grid(x, y))
 op <- par(mfrow = c(2,1), mar = .1 + c(3,3,0,0)) # グラフィックス画面を上下に二分割
 contour(x, y, matrix(log10(z), length(x)))       # 曲面の等高線を眺める
 nlm.f2 <- nlm(fR, c(-1.2, 1), hessian=TRUE)      # 最小化実行
 points(rbind(nlm.f2$estim[1:2]), col = "red", pch = 20)  # 等高線図に最適パラメータをプロット
 ## 拡大曲面を画く
 rect(.9, .9, 1.1, 1.1, border = "orange", lwd = 2)
 x <- y <- seq(.9, 1.1, len=100)
 z <- fx(expand.grid(x, y))
 contour(x, y, matrix(log10(z), length(x)))   # 曲面の等高線を眺める
 mtext("zoomed in")                           # 注釈を入れる
 box(col = "orange")                          # オレンジ色の枠を描
 points(rbind(nlm.f2$estim[1:2]), col = "red", pch = 20)   # 等高線図に最適パラメータをプロット
 par(op)       # グラフィックスパラメータ復帰
 nlm.f2
}
> nlm.f2 <- nlm2()
> nlm.f2
$minimum           # 真値は当然ゼロ
[1] 3.973766e-12
$estimate          # 真のパラメータは当然 c(1,1) 
[1] 0.999998 0.999996
$gradient
[1] -6.539279e-07  3.335998e-07
$hessian
          [,1]      [,2]
[1,]  802.2368 -400.0192
[2,] -400.0192  200.0000
$code  # 満足のいく終了
[1] 1
$iterations
[1] 23
> str(nlm.f2)
List of 6
 $ minimum   : num 3.97e-12
 $ estimate  : num [1:2] 1 1
 $ gradient  : num [1:2] -6.54e-07  3.34e-07
 $ hessian   : num [1:2, 1:2]  802 -400 -400  200
 $ code      : int 1
 $ iterations: int 23
> det(nlm.f2$hessian) # hessian の行列式値(正、つまり正定値符合行列!)
[1] 432.0022 

#ref(): The style ref(filename,pagename) is ambiguous and become obsolete. Please try ref(pagename/filename)

例3 demo(nlm) より (gradient を与える例)

nlm3 <- function () {
  fg <- function(x) {  # x は二次元ベクトル変数 (x[1],x[2])
    gr <- function(x1, x2) { # fg の偏微分(gradient)ベクトル関数
        c(-400*x1*(x2 - x1*x1)-2*(1-x1), 200*(x2 - x1*x1))
    }
    x1 <- x[1]; x2 <- x[2]
    res<- 100*(x2 - x1*x1)^2 + (1-x1)^2 # fg の本体
    attr(res, "gradient") <- gr(x1, x2) # fg の gradient 属性(関数本体)は gr であることを教える
    return(res)
  }
  nlm(fg, c(-1.2,1), hessian=TRUE) # 非線形回帰実行
}
> nlm3.f <- nlm3() # 非線形回帰実行
> nlm3.f           # 結果表示
$minimum           # 最小値
[1] 1.182096e-20
$estimate          # その時のパラメータ値 
[1] 1 1 
$gradient          # その時の gradient 値 (数値誤差範囲内で零に近い!)
[1]  2.583521e-09 -1.201128e-09
$hessian           # その時の hessian 値
        [,1]    [,2]
[1,]  802.24 -400.02
[2,] -400.02  200.00
$code              # ベストな終了判定コード
[1] 1
$iterations        # 収束までに24回の繰返しが行なわれた
[1] 24

> str(nlm3.f)  # 簡潔な要約
List of 6
 $ minimum   : num 1.18e-20
 $ estimate  : num [1:2] 1 1
 $ gradient  : num [1:2]  2.58e-09 -1.20e-09
 $ hessian   : num [1:2, 1:2]  802 -400 -400  200
 $ code      : int 1
 $ iterations: int 24
# 目的関数の内容
# gradient の関数本体が属性として fg の定義の一部を構成することを注意。
> fg          
function(x) {  # x は二次元ベクトル変数 (x[1],x[2])
    gr <- function(x1, x2) { # fg の偏微分(gradient)ベクトル関数
        c(-400*x1*(x2 - x1*x1)-2*(1-x1), 200*(x2 - x1*x1))
    }
    x1 <- x[1]; x2 <- x[2]
    res<- 100*(x2 - x1*x1)^2 + (1-x1)^2 # fg の本体
    attr(res, "gradient") <- gr(x1, x2) 
    return(res)
 }

例4 (例3と同じ、gradient を数式微分で計算する例)

fdd は gradient 関数そのものではなく、gradient 属性を備えた目的関数自身であることを注意。

> nlm4 <- function () { 
    fd <- deriv(~ 100*(x2 - x1*x1)^2 + (1-x1)^2, c("x1", "x2")) #  目的関数の数式微分
    fdd <- function(x1, x2) {}  # 空の関数定義
    body(fdd) <- fd             # 関数の本体を fd とする
    nlm(function(x) fdd(x[1], x[2]), c(-1.2,1), hessian=TRUE) 
}
> nlm4.f <- nlm4()
> nlm4.f
$minimum
[1] 1.182096e-20
$estimate
[1] 1 1
$gradient
[1]  2.583521e-09 -1.201128e-09
$hessian
        [,1]    [,2]
[1,]  802.24 -400.02
[2,] -400.02  200.00
$code
[1] 1
$iterations
[1] 24
> str(nlm4.f)
List of 6
 $ minimum   : num 1.18e-20
 $ estimate  : num [1:2] 1 1
 $ gradient  : num [1:2]  2.58e-09 -1.20e-09
 $ hessian   : num [1:2, 1:2]  802 -400 -400  200
 $ code      : int 1
 $ iterations: int 24
# fdd は gradient 関数そのものではなく、gradient 属性を備えた目的関数自身であることを注意。
# (オブジェクト名の先頭のドットは紛らわしいオブジェクト名を避けるため)
> fdd
function (x1, x2)
{
    .expr2 <- x2 - x1 * x1
    .expr5 <- 1 - x1
    .value <- 100 * .expr2^2 + .expr5^2
    .grad <- array(0, c(length(.value), 2), list(NULL, c("x1",
        "x2")))
    .grad[, "x1"] <- -(2 * .expr5 + 100 * (2 * ((x1 + x1) * .expr2)))
    .grad[, "x2"] <- 100 * (2 * .expr2)
    attr(.value, "gradient") <- .grad
    .value
}

例5 demo(nlm) より (gradient と hessian をともに式で与える(失敗する)例)

hessian を陽に与えてもそれだけ賢くなるわけではない!? 数値微分の値は正確な値とは当然微妙にことなるはずなので、こういうこともおき得る!?

これはhessianの定義が不完全(もしくはmatrixでの行と列の入れ間違い)だからでしょうか? 単なる誤植だったとすると、微分を明確に与えることによって数値微分よりも 曲面が複雑になってしまったからでしょうか?

hessian の定義は正しいですね(どのみち対称ですから)。失敗する理由は良くわかりませんが、真の解 (1,1) における hessian 行列の固有値は 1001.6 と 0.39936 と極めてアンバランス(ちょうどデスバレー風)なのが原因では?同じ変位でも方向により関数値が極端に変わります。

> fgh <- function(x)
  {
    gr <- function(x1, x2) {
        c(-400*x1*(x2 - x1*x1) - 2*(1-x1), 200*(x2 - x1*x1))
    }
    h <- function(x1, x2) {
        a11 <- 2 - 400*(x2 - x1*x1) + 800*x1*x1
        a21 <- -400*x1
        matrix(c(a11,a21,a21,200),2,2)
    }
    x1 <- x[1]; x2 <- x[2]
    res<- 100*(x2 - x1*x1)^2 + (1-x1)^2
    attr(res, "gradient") <- gr(x1, x2) # fgh に gradient 属性として gr を与える
    attr(res, "hessian") <- h(x1, x2)   # fgh に hessian 属性として h を与える
    return(res)
}
> nlm5.f <- nlm(fgh, c(-1.2,1), hessian=TRUE)
> nlm5.f
$minimum        # 完全な失敗!
[1] 2.829175
$estimate
[1] -0.6786981  0.4711891
$gradient
[1] -0.4911201  2.1115987
$hessian
         [,1]     [,2]
[1,] 366.1188 271.4593
[2,] 271.4593 200.0000
$code           # 失敗!
[1] 4
$iterations     # 最大繰返し回数まで実行されている! 
[1] 100
> eigen(nlm5.f$hessian)  # hessian は正定値符合行列でない!
$values
[1] 566.9414302  -0.8225813
$vectors
          [,1]       [,2]
[1,] 0.8039230 -0.5947334
[2,] 0.5947334  0.8039230
> fgh               # 関数 fgh の中身は
function(x)
{
    gr <- function(x1, x2) {
        c(-400*x1*(x2 - x1*x1) - 2*(1-x1), 200*(x2 - x1*x1))
    }
    h <- function(x1, x2) {
        a11 <- 2 - 400*(x2 - x1*x1) + 800*x1*x1
        a21 <- -400*x1
        matrix(c(a11,a21,a21,200),2,2)
    }
    x1 <- x[1]; x2 <- x[2]
    res<- 100*(x2 - x1*x1)^2 + (1-x1)^2
    attr(res, "gradient") <- gr(x1, x2)
    attr(res, "hessian") <- h(x1, x2)
    return(res)
}
# 反復回数を増やせばそれなりの結果が得られる
> nlm(fgh, c(-1.2,1), iterlim=10000, hessian=TRUE)
$minimum
[1] 1.675511e-07
$estimate
[1] 0.9995911 0.9991805
$gradient
[1] -0.0000996657 -0.0003592498
$hessian
          [,1]      [,2]
[1,]  801.5865 -399.8564
[2,] -399.8564  200.0000
$code
[1] 2
$iterations
[1] 5044

例6 数式微分関数 deriv3 でgradient と hessian を同時に与える

> fd <- deriv3(~ 100*(x2 - x1*x1)^2 + (1-x1)^2, c("x1", "x2")) #  関数の数式微分
> fdd <- function(x1, x2) {}  # 空の関数定義
> body(fdd) <- fd             # 関数の本体を fd とする

> fdd   # fdd の中身 (gradient と hessian がともに属性として備わっている)
function (x1, x2)
{
    .expr2 <- x2 - x1 * x1
    .expr5 <- 1 - x1
    .expr9 <- x1 + x1
    .value <- 100 * .expr2^2 + .expr5^2
    .grad <- array(0, c(length(.value), 2), list(NULL, c("x1",
        "x2")))
    .hessian <- array(0, c(length(.value), 2, 2), list(NULL,
        c("x1", "x2"), c("x1", "x2")))
    .grad[, "x1"] <- -(2 * .expr5 + 100 * (2 * (.expr9 * .expr2)))
    .hessian[, "x1", "x1"] <- -(100 * (2 * ((1 + 1) * .expr2 -
        .expr9 * .expr9)) - 2)
    .hessian[, "x1", "x2"] <- .hessian[, "x2", "x1"] <- -(100 *
        (2 * .expr9))
    .grad[, "x2"] <- 100 * (2 * .expr2)
    .hessian[, "x2", "x2"] <- 100 * 2
    attr(.value, "gradient") <- .grad
    attr(.value, "hessian") <- .hessian
    .value
}

例7 example(nlm) より (追加引数つきの目的関数の例)

> f <- function(x, a) {  # ベクトル化関数 (変数次元は未設定、a はパラメータ)
   res <- sum((x - a)^2)
   attr(res, "gradient") <- 2 * (x - a)
   res
}
# 目的関数は二変数関数 (x[1]-3)^2+(x[2]-5)^2 とされる。
> nlm(f, c(10, 10), a = c(3, 5)) # 二次元引数 x、a はパラメータ
$minimum  # 当然
[1] 0
$estimate # =a、当然
[1] 3 5
$gradient
[1] 0 0
$code
[1] 1
$iterations 
[1] 1

添付ファイル: filenlm1.png 1784件 [詳細] filenlm2.png 1860件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2023-03-25 (土) 11:19:17