1

Size: px

Start display at page:

Download "1"

ときよせ
5 years ago
Views:

1 R による非線形最小二乗法. 非線形回帰モデル回帰モデルにおいて被説明変数が未知パラメータについて線形である場合は線形回帰モデルとよばれる例えば以下のようなモデルはすべて線形回帰モデルの例である ( 例 ) y x, ( 例 2) y log( x ), ( 例 3) log y x このようなモデルの場合は通常の最小二乗法によって未知パラメータ α β を推定する事ができるこのような線形回帰モデルに対して被説明変数が未知パラメータについて非線形である場合は非線形回帰モデルと呼ばれる例えば以下のようなモデルである ( 例 4) y exp( x ), ( 例 5) y exp( x ) ( 例 6) y sn( x ), ( 例 7) y x z このようなモデルはより一般的に y m( x, θ) と表せるここで x は説明変数のベクトル θ は未知パラメータのベクトル ε は直接観測できない誤差項である説明変数 x と誤差項 ε は独立であると仮定する例えば上記の例 5 では x = x, θ = (α, β), m( x, θ ) exp( x ) となるこのような非線形回帰モデルに含まれる未知パラメータを推定するには非線形最小二乗法と呼ばれる方法を用いる非線形最小二乗法の具体的なやり方はここでは述べない ( 線形の最小二乗法と原理的にはまったく同じである ) 以下では R によって上記のようなモデルを推定する方法のみを述べる 2. R 関数 nls による非線形最小二乗法 2. 人口データの非線形モデル推定 R によって非線形最小二乗推定法を行うやり方を簡単に説明する説明のために用いるデータとしてパッケージの car の中に入っているアメリカの 0 年ごとの人口の推移のデータ ( 単位 : 百万人 ) を使用するまずパッケージ car をダウンロード ( インスツール ) する ( パッケージのダウンロードおよびインスツールの仕方はすでに説明したので省略 ) 次に > lbrary(car) としてパッケージを R に読み込ませて使用できるようにするパッケージを読み込んだのでここに入っている USPop というデータ ( アメリカの人口推移 ) が使用できる最初の 5 行を見てみると > head(uspop,5) populaton となっているは西暦であるデータの特徴を視覚的にとらえるために散布図を描いてみる > plot(populaton~, data=uspop, man="u.s. populaton") この資料は私のゼミおよび講義で R の使用法を説明するために作成した資料ですホームページ上で公開しており自由に参照して頂いて構いませんただし内容について一応検証してありますが間違いがあるかもしれません間違いがあった場合でもそれによって生じるいかなる損害不利益について責任は負いかねますのでご了承ください

2 すると以下のような散布図が描かれる U.S. populaton populaton この図より明らかなように人口は西暦に比例して増えていくもののその関係は直線的ではない例えばここに最小二乗法によって推定した直線を書き込んでみよう > ablne(lm(populaton~,data=uspop)) すると以下のように直線が書き込まれるこの直線はあまりデータにフィットしていないことがわかる U.S. populaton populaton 上記から考えられるのは人口の推移は西暦と非線形の関係にあるであろうという事であるそこで以下のような非線形回帰モデルを考える y m( x, ) exp[ ( 0, =,,N. x )] 2

3 ここで y は人口 x は西暦である未知パラメータは θ = (α, β 0, β ) の 3 つであるこのような非線形のモデルを推定する R の関数として nls 関数がある上記のモデルを nls 関数を用いて推定してみよう以下のように入力する > result=nls(populaton~a/(+exp(-(b0+b*))), + start=lst(a=400,b0=-49,b=0.025),data= USPop, trace=true) ここで populaton は USPop のデータにある populaton が被説明変数である事を表している ~ の後の a/(+exp(-(b0+b*))) はそれに対して具体的にどのような非線形の関数をあてはめるかを指定しているは USPop にあるという変数であるコマンドを打ち込んでいる途中で Enter キーを押せば次の行に移りその際 + と表示されるここでは上記の部分を打ち込んだ後 Enter キーによって次の行に移っているそこから次の start=lst(a=400,b0=-49,b=0.025) を打ち込んでいるがこれは初期値と呼ばれるものである非線形最小二乗法を始めるにはまず未知パラメータの値として初期値を指定する必要があるそこから収束計算によって未知パラメータの真の値の推定値を計算する初期値は基本的には適当な値でよいが通常この初期値が真の値に近ければ近いほど収束計算が早く終わる逆に初期値が真の値から離れているほど計算に時間がかかるまた場合によっては初期値が異なると計算結果が異なるという事もあるいくつか異なった初期値を試すのが普通であるここでは上記のようにした最後の trace=true は途中の収束計算の経過を表示するかどうかを指定するものである trace=true であれば途中の収束計算の経過を表示し trace=false であれば表示しない何も入力しなければ自動的に trace=false となる Enter キーを押して上記のコマンドを実行すると計算の途中結果が表示される左からそれぞれ非線形最小二乗法の残差平方和の値 a の値 b0 の値 b の値の収束計算の途中経過が表示される推定結果は summary 関数によって見ることができる > summary(result) Formula: populaton ~ a/( + exp(-(b0 + b * ))) Parameters: Estmate Std. Error t value Pr(> t ) a e-0 *** b e-5 *** b e-5 *** --- Sgnf. codes: 0 *** 0.00 ** 0.0 * Resdual standard error: on 9 degrees of freedom Number of teratons to convergence: 6 Acheved convergence tolerance:.48e-06 パラメータの推定値は Parameters: の下の Estmate の列に並んでいる Std.Error, t value, Pr(> t ) は標準誤差それぞれ t 値 ( 推定値 / 標準誤差 : 真の値が 0 という帰無仮説の検定用 ) その P 値である Resdual standard error は誤差項 ε の標準偏差の推定値である推定値の下でのこの非線形回帰モデルの残差は > resduals(result) によって見ることができる ( 確認してください ) また推定値の下でのあてはめ値は 3

4 > predct(result) によって見ることができる ( 確認してください ) 実際の西暦と人口の値の散布図に非線形の当てはめ値の曲線を書き込むこともできる ( 以前の図を閉じて ) 改めて実際の西暦と人口の散布図を書いておこう > plot(populaton~,uspop) 次にこの散布図に当てはま値の曲線を加えるために > lnes(seq(790,2000,by=0),predct(result)) と入力しエンターキーを押すここで seq(790,2000,by=0) は 790 から 2000 まで 0 ずつ増加させた数列 ( 英語で sequence) のデータを作るという意味であるすると散布図に当てはめ値の曲線が書き込まれた以下の図が出力される populaton 上記の手順は例えば曲線を書き込む手順をいくつかに段階わけして > pre=predct(result) > p=seq(790,2000,by=0) ( ここで p は予測 -predct- 用の西暦という事で p を付けているが名前は何でもよい ) > lnes(p,pre) としても同じ図が出力される ( 確認してください ) さらに残差をプロットするには > res=resduals(result) > plot(uspop$,res,type='b') > ablne(h=0,lty=2) というコマンドを実行する 2 行目の ablne は線を加えるという事で h=0 は y 軸の 0 の部分に水平な (horzontal) 線を加えるという事である lty=2 はその線が点線となるようにするという事である以下のような図が出力される 4

5 USPop$ またもとのデータは 2000 年までしかないがそれ以降の年の人口の予測値を推定値の下でのこの非線形モデルを使用して予測したものを計算することもできるまず > p2 = data.frame(=seq(790,200,by=0)) によって 790 年から 200 年までの 0 年ごとのデータを作る次に > pre2=predct(result,p2) とすれば 200 年以降 (23 個めから ) の予測値も出力されそれに pre2 という名前がつけられる最後に予測図と実際の値との図は > plot(populaton~,uspop,xlm=c(790,200),ylm=c(0,450)) > lnes(p2$,pre2) によって描くことができる以下のような図が出力される populaton res

6 2000 年以降の予測値が図示されているのがわかる 2.2 誤差項を AR() モデルにしてみる上の残差のプロットは残差に正の系列相関があることを示唆している実際に残差の自己相関を > acf(resduals(result),lag.max=0) で計算してみると以下のようになる Seres resduals(result) ACF これより次の自己相関が高いことがわかるよって ( この図ではむしろ MA() のように見えるがひとます ) 誤差項が AR() モデルに従っているとしてモデル化してみよう以下のようなモデルを考える y exp[ ( 0 x )], u, u ~..d. (0, σ 2 ) ここで ε = y α /(+exp[ (β 0 + β x )]) を ε = ε + u に代入するとこのモデルは y y u () exp[ ( 0 x )] exp[ ( 0 x )] と表すことができるこのモデルを推定してみようこのモデルの推定において注意すべき点のつがの推定である定常性を仮定するならの値の取りうる範囲は < < となるよって推定値もこの範囲の値にならなければならないが通常統計ソフトの最適化アルゴリズムでは未知パラメーターについて最適化する際パラメーターの取りうる範囲に制約を課さずに最適化するものがほとんどであるつまりこの場合そのようなアルゴリズムを用いて上記のを推定すると計算の途中でが上記の範囲の外側の値を取ってしまう可能性がありその場合は目的関数がおかしな値をとりアルゴリズムがストップしてしまう ( ただし後に見るように実は nls 関数にはパラメーターの範囲を指定することができるアルゴリズムが存在するがひとまずそれは使用しないものとして話を進める ) このような問題を避けるためには例えばを直接推定するのではなくはある無制約パラメーター κ の関数でその関数はからの間の値をとるすなわち ( ) ( ) であるして κ を推定し κ の推定値 ˆ を用いてを ˆ ( ˆ ) と推定する方法が考えられるこのような変換としては例えば ( ) 2( ), Lag ( ) exp( ) が考えられる ( ) はロジスティック分布の分布関数であり (0, ) の範囲の値を取るがそれを 2 倍してを引いているので ( ) は (, ) の範囲の値を取ることになるそれでは () 式のモデルを推定するただし今回は自己相関は正の値なのでの取りうる値としては (0, ) とすれば十分であろうよって ( ) ( ) とするこの時 () 式をパラメーター κ について書き直 6

7 すと y y u (2) exp[ ( 0 x )] exp( ) exp[ ( 0 x )] となるこのモデルの推定において必要な系列は y y, x, x であるのでそれらを作っておこう ( データの数が一つ減ることに注意 ) 元のデータを見てみると > USPop$populaton [] [6] [] [6] [2] > USPop$ [] [3] なので元のデータはそれぞれ 22 個あるここからまず y の系列を作る ( これを pop とする ) > pop=uspop$populaton[2:22] > pop [] [6] [] [6] [2] 同様に y の系列は ( これを popl とする L は lag の L) > popl=uspop$populaton[:2] > popl [] [6] [] [6] [2] 同様に x と x の系列は > =USPop$[2:22] > L=USPop$[:2] と作ることができるそれでは (2) 式を推定してみよう以下のコマンドを打ち込む ( 少し長いが ) > result2=nls(pop~a/(+exp(-(b0+b*))) + +(/(+exp(-k)))*(popl-(a/(+exp(-(b0+b*l))))), + start=lst(a=400,b0=-49,b=0.025,k=0),trace=true) ( 適当なところで改行していることに注意 ) 結果は > summary(result2) Formula: pop ~ a/( + exp(-(b0 + b * ))) + (/( + exp(-k))) * (popl - (a/( + exp(-(b0 + b * L))))) Parameters: Estmate Std. Error t value Pr(> t ) a * b e-05 *** 7

8 b *** k Sgnf. codes: 0 *** 0.00 ** 0.0 * Resdual standard error: 3.42 on 7 degrees of freedom Number of teratons to convergence: 2 Acheved convergence tolerance: 6.28e-06 となる ( ここで κ の推定値が有意になっていないがこれは直観的にはが 0.5 と有意に異ならないということでが 0 と異ならないということではないことに注意 ) κ の推定値からの推定値を計算してみると > /(+exp( )) [] であるだいたいˆ 0.84 でありかなり高い値になっている残差をプロットしてみると > plot(resduals(result2),type="l") resduals(result2) Index のようになるその自己相関は > acf(resduals(result2),lag.max=0) Seres resduals(result2) ACF となる一応自己相関が消えているのがわかるまたあてはめ値のグラフは > plot(pop~) > lnes(seq(800,2000,by=0),predct(result2)) Lag より下の図のようになる (pop とは 800 から 2000 までの値であることに注意 ) 下の図のさらに 8

9 下の図はもとのあてはめ値のグラフ populaton pop 若干あてはまりがよくなっているのが見て取れるこれは推定結果のところの誤差項の分散が小さくなっていることからもわかるこのモデルを用いて予測を行うこともできるが若干の注意が必要であるこのモデルでは説明変数に被説明変数のラグが入っているため期先予測であれば ( 被説明変数の ) 現在の値 2 期先予測であれば期先の値が必要になる期先の値はもちろんまだ観測されていないのでその予測値で置き換えて 2 期先予測を行いまたその 2 期先予測値を用いて 3 期先予測を行うという逐次的に予測値を計算する必要が出てくる 4. アルゴリズム port を使用してパラメーターの範囲を指定するさて上記では最適化のアルゴリズムにおいてパラメーターの範囲が指定できない場合の対処法を見てみたが実は nls 関数では最適化の際にパラメーターの範囲を指定できる port というアルゴリズムが使用できるこのアルゴリズムを使用することによって最適化の際にそれぞれのパラメーターが取りうる範囲の上限と下限を指定することができるこのアルゴリズムを用いて () 式のモデルにおいてを直接推定してみよう ( ここでは f というパラメーターで表すとする ) 以下のように入力する > result2=nls(pop~a/(+exp(-(b0+b*)))+ + f*(popl-(a/(+exp(-(b0+b*l))))), + alg="port",start=lst(a=400,b0=-49,b=0.025,f=0), + lower=lst(a=0,b0=-inf,b=0,f=-),upper=lst(a=inf,b0=inf,b=inf,f=)) 9

10 ここで alg="port" によって port ( 何かの略だろうがよくわからない ) というアルゴリズムを使用することを指定し lower=lst(a=0,b0=-inf,b=0,f=-) でそれぞれのパラメーターの下限を upper=lst(a=inf,b0=inf,b=inf,f=)) で上限を指定している (Inf はを意味する ) 推定結果は > summary(result2) Formula: pop ~ a/( + exp(-(b0+b*)))+f*(popl-(a/(+exp(-(b0 + b * L))))) Parameters: Estmate Std. Error t value Pr(> t ) a * b e-05 *** b *** f *** --- Sgnf. codes: 0 *** 0.00 ** 0.0 * Resdual standard error: 3.42 on 7 degrees of freedom Algorthm "port", convergence message: relatve convergence (4) となる今回も f の推定値がほぼ 0.84 でありそのほかの推定値もほぼ同じである 0

> usdata01 と打ち込んでエンターキーを押すと V1 V2 V : : : : のように表示され読み込まれていることがわかるここで V1, V2, V3 は R が列のデータに自動的につけた変数名である ( variable

> usdata01 と打ち込んでエンターキーを押すと V1 V2 V : : : : のように表示され読み込まれていることがわかるここで V1, V2, V3 は R が列のデータに自動的につけた変数名である ( variable R による回帰分析 ( 最小二乗法 ) この資料では 1. データを読み込む 2. 最小二乗法によってパラメーターを推定する 3. データをプロットし回帰直線を書き込む 4. いろいろなデータの読み込み方について簡単に説明する 1. データを読み込む以下では read.table( ) 関数を使ってテキストファイル ( 拡張子が.txt のファイル ) のデータの読み込み方を説明する 1.1