Levenberg-Marquardtアルゴリズムで非線形回帰分析
Rでminpack.lmパッケージのnls.lm関数を使うとLevenberg-Marquardt法で非線形回帰分析を行うことができる。言いかたはマルカート法、マーカート法？、フランスの方のようなのでマーカール法とも言われているようだ。私はマルカートで覚えていた。Levenberg-Marquardt法は非線形最小二乗問題を解く手法として広く使われている。最急降下法・ニュートン法を組み合わせた方法で現在の解が正解から遠い場合は遅いが収束することが保証されている最急降下法と同じように動作し、正解から近い場合はニュートン法を実行するとのこと。
非線形回帰を行うことになった背景
最近ECサイトの分析を行うことが多く、よく言われているベキ乗則を随所に見かける。データ例
> head(d) X Y 1 2 9794688 2 3 5973376 3 9 2946944 4 10 2649600 5 12 2520064 6 15 2168766 > plot(d$X, d$Y)
これは、とあるカテゴリーのランキングと指標○○の関係。こんな関数で非線形回帰分析を行う。
yは指標○○、xはランキング、aは切片(1位の指標○○値)、bは減衰項
nls.lm関数を実行する
> library(minpack.lm) > res_lm <- lm(formula=log(Y)~log(X), data=d) # 初期パラメータ設定の為に対数変換して線形回帰 > lm_pa <- exp(res_lm$coefficients[1]) # aの初期パラメータ > lm_pb <- res_lm$coefficients[2] # bの初期パラメータ > x <- d$X # 説明変数 > obs <- d$Y # 目的変数 > pred <- function(parS, xx) parS$a * xx ^ parS$b # モデル式 > resid <- function(p, observed, xx) observed - pred(p,xx) # 残差 > parStart <- list(a=lm_pa, b=lm_pb) # 初期パラメータ > nls.out <- nls.lm(par=parStart, fn=resid, observed=obs, xx=x, control=nls.lm.control(maxiter=1024,nprint=1)) It. 0, RSS = 1.59819e+15, Par. = 1.29507e+08 -1.51893 It. 1, RSS = 1.2133e+14, Par. = 3.83967e+06 -1.26897 It. 2, RSS = 9.01298e+12, Par. = 2.451e+07 -1.18445 It. 3, RSS = 4.13613e+12, Par. = 1.32097e+07 -0.719333 It. 4, RSS = 1.38756e+12, Par. = 1.64626e+07 -0.823374 It. 5, RSS = 1.34457e+12, Par. = 1.66165e+07 -0.817677 It. 6, RSS = 1.34456e+12, Par. = 1.66199e+07 -0.817936 It. 7, RSS = 1.34456e+12, Par. = 1.66198e+07 -0.81793 > nls.out Nonlinear regression via the Levenberg-Marquardt algorithm parameter estimates: 16619797.5617096, -0.817929785073928 residual sum-of-squares: 1.345e+12 reason terminated: Relative error in the sum of squares is at most `ftol'. > > nls_a <- nls.out$par$a > nls_b <- nls.out$par$b > cor(nls_a*x^nls_b, obs) [1] 0.995052 > > plot(d$X, d$Y) > x <- seq(min(d$X), max(d$X), by=1) > x.pred <- nls_a*x^nls_b > lines(x, x.pred, col=2) >
あてはまり度合も良好で説明力も高そう。今回の非線形回帰の初期値はモデル式を対数変換して線形回帰した回帰係数を使用した。Levenberg-Marquardt法で得られた回帰係数はRで一般的に知られているnls関数で得られる値とほとんど変わらない。
nlsだとエラーが出る場合が。。。
Levenberg-Marquardt法を行う前はnls関数を使用していた。nls関数の場合は分析するデータにより下記のようにエラーが出る場合がある。algorithm="port"で改善する場合はあったが全てのエラーは解消しなかった。今のところnls.lm関数ではエラーは出てない。
Error in nls(formula = Y ~ a * X^b, data = d, start = list(a = exp(res_lm$coefficients[1]), : 勾配が特異です Error in numericDeriv(form[[3L]], names(ind), env) : モデルを評価する際，欠損値または無限大が生成されました
参考
