読者です 読者をやめる 読者になる 読者になる

でたぁっ 感動と失敗の備忘録

データ解析を担当することになったサラリーマンの備忘録

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)

f:id:deta:20140508043056p:plain
これは、とあるカテゴリーのランキングと指標○○の関係。こんな関数で非線形回帰分析を行う。
y=ax^b
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)
> 

f:id:deta:20140508044224p:plain
あてはまり度合も良好で説明力も高そう。今回の非線形回帰の初期値はモデル式を対数変換して線形回帰した回帰係数を使用した。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) : 
   モデルを評価する際,欠損値または無限大が生成されました 

参考