ワイブル分布を使った推定

質問

ある製品が販売後t年経っても廃棄されずに残っている率（残存率Y(t)と呼ぶことにします）のデータがt=1からt=10まであります。これを元に、t=11からt=50までの残存率を推定したいと思っています。

ワイブル分布を説明した
http://ja.wikipedia.org/wiki/%E3%83%AF%E3%82%A4%E3%83%96%E3%83%AB%E5%88%86%E5%B8%83
を見て、ワイブル・プロットという方法を使うと次のようにできると考えたのですが、妥当なやり方でしょうか。

・t=1～10の時のln(t)を説明変数、ln(ln(1/Y(t)))を被説明変数とした単回帰を行い、回帰式を求める
・回帰式にt=11～50を代入して、ln(ln(1/Y(t)))の推定値を求める（これをZ(t)と書くことにする）
・1/(exp(exp(Z(t)))を計算することにより、Y(t)を求める。

Excelで作ってみたところ、それっぽい曲線はできたのですが。

よろしくお願いいたします。

回答No.3

stomachman · ベストアンサー

ANo.2へのコメントについて。

> 線形回帰の関数slope, intercept等

　実はそんなの使ったことなくてですね、大抵いつも行列計算の関数でやってるんですよ。

　要するに
　　r[n] = a x[n] + b - y[n] 　　(n=1,2,…,N。a, bが未知数、r[n]は残差、x[n],y[n]は値が分かっている。)
においてΣ((w[n]r[n])^2) が最小になる解 (a,b)を計算したい。（ただしw[n]は値が分かっている。）ということは
　　 w[n]r[n] = a w[n]x[n] + w[n]b - w[n]y[n]　　(n=1,2,…,N)
においてΣ((w[n]r[n])^2) が最小になる解 (a,b)を計算したって同じなわけで、ここで
　　Y[n] = w[n]y[n]
　　X[n] = w[n]x[n]
　　R[n] = w[n]r[n]
と書くことにすれば
　　R[n] = a X[n] + b w[n] - Y[n]　　(n=1,2,…,N)
においてΣ(R[n]^2) が最小になる解 (a,b)を計算すればよし。

Excelでやるんですと、

A1:ANにX[n]を入力。B1:BNにw[n]を入力。　つまり、一次式の係数を並べたN行2列の行列（これをJacobianと言いますんで、Jとしましょう）をA1:BNに作る。

C1:CNにY[n]を入力。

E1:G2を選択した状態でE1に {=mmult(transpose(A1:BN), A1:CN)}　これでE1:F2に J' Jを、G1:G2にはJ' Y を入れたことになります。（ただし、 ' は行列の転置という意味です。）

I1:J2を選択した状態でI1に {=minverse(E1:F2)}　これでI1:J2に(J' J)の逆行列ができます。

L1:L2を選択した状態でL1に{=mmult(I1:J2,G1:G2)}　すなわち、(J' J)の逆行列とJ' Yとの積がL1:L2に入ります。

で、答はL1=a, L2=b でやんす。

　なお、{ } ってのは、その内側の式をセルに入力しといて、数式の入力欄内にカーソルを入れた状態のまんまで [command]キーと[enter]キーを一緒に押す（のがMac。Windowsだとどうやるんだか忘れた。）と、勝手に{ } が付きます。

　この手なら、パラメータがもっと沢山ある場合の線形回帰もできますから、知ってて損はないでしょう。

お礼

ありがとうございます！
できました。
元の方法よりも、残差平方和が35%ほど小さくなりました。

詳しくお教えいただきまして、本当にありがとうございました。

回答No.4

stomachman

ついでに非線形最小二乗法の数値計算のやりかたも書いときましょう。大変そうに思うでしょうけど、実はそうでもないんです。

　タテ軸のデータy[n]に対して、パラメータp,qとヨコ軸の値t[n]で値が決まるような既知の非線形関数f(p,q,t[n])による近似の残差が
　　r[n] = f(p,q,t[n]) - y[n] (n=1,2, …, N）
であり、しかも、Σr[n]^2を（最小とは言えないまでも）かなり小さくするp, qの値p0, q0（つまり、最小二乗解の近似値）が分かっているとします。
　　r0[n] = f(p0,q0,t[n]) - y[n] (n=1,2, …, N）
で、これを改良した、もっと良い近似解p1=p0+Δp, q1=p0+Δqが欲しい。

　そこで、fをpで偏微分して得られる関数を(∂f/∂p)(p,q,t)、(∂f/∂q)(p,q,t)も同様としまして、1次までのテイラー展開でf(p1,q1,t)を近似すると、
　　f(p1,q1,t) = f(p0+Δp,q0+Δq,t) ≒ f(p0,q0,t) + ((∂f/∂p)(p0,q0,t))Δp + ((∂f/∂q)(p0,q0,t))Δq
この式は、各tについて、(p,q,f(p,q,t))のグラフ（曲面）を描いた時の点(p0,q0)におけるf(p,q,t)の接平面を表しています。これを使うと
　　r[n] ≒ ((∂f/∂p)(p0,q0,t[n]))Δp + ((∂f/∂q)(p0,q0,t[n]))Δq + r0[n]
となり、右辺は既知の数値((∂f/∂p)(p0,q0,t[n])), ((∂f/∂q)(p0,q0,t[n]))を係数とし、既知の定数r0[n]を定数項とした、未知数Δp, Δqの一次式になっています。ですから、この右辺を使ってΣr[n]^2を(近似的に)最小にするΔp, Δqを計算すれば、最小二乗解p0,q0を改良できる。

　手順は結局、Jacobianとして
A1:ANに((∂f/∂p)(p0,q0,t[n]))を入力。
B1:BNに((∂f/∂q)(p0,q0,t[n]))を入力。　
そして、
C1:CNに-r0[n]を入力。
　あとはANo.3と同じ操作で（a, bの代わりに）Δp, Δqが得られますから、p0をp1=p0+Δpで、q0をq1=q0+Δqでそれぞれ置き換えて改良された解を得る。
　これを繰り返せば、最小二乗解に近づいて行ける訳です。

お礼

ありがとうございます。
繰り返しを手作業で行うのであれば、おっしゃるとおり操作はそんなに大変ではないですね。

回答No.2

stomachman

ANo.1へのコメント、了解です。

　fittingと回帰は同義でしょう。
　問題は、線形回帰で計算することの適切さであり、何をやりたいかによって、そのやり方が適切かどうかが決まります。何をズレの度合い（残差）と見て、その残差をどういうやりかたで評価するのが適切か、ということです。
　まずは残差について、「最終的に見たい尺度で描いたグラフ」（たとえばご質問にあるグラフだとか、縦軸がlogのグラフ）におけるYとZの差（タテ方向のずれ）を残差ε(t) と考えるのが適切でしょう。また残差の評価関数は、「全体が合う感じにしたい」という程度のことであれば、二乗和E=Σ(ε(t))^2 を使うのが簡単です。
　だとすれば、
(a) ご質問にあるグラフなら残差
　　ε(t) = Y(t)-Z(t)
の二乗和
　　E=Σ(ε(t))^2　　(Σはt=t1,t2,…,tnの総和)
が最小になるようにZのパラメータを決めるのが良いし、
(b) 縦軸が対数のグラフ（つまり残ってる個数の比率を見たい）なら残差
　　δ(t) = ln(Y(t))-ln(Z(t))
の二乗和
　　D=Σ(δ(t))^2
の最小化を考えるのが良いでしょう。

　ところが、
(c)ご質問の計算では、残差
　　γ(t) = ln(-ln(Y(t)))-ln(-ln(Z(t)))
の二乗和
　　G=Σ(γ(t))^2
が最小になるようにZのパラメータを決めています。なので、Yが1に近いときと、Yが0に近いときにln(-ln(Z(t))はそのデータは重視される。なぜなら、Zをわずかに変えただけでln(-ln(Z(t)))が大きく変化するから、変化がγ(t)に大きく反映する訳です。またY=1/eに近いとき、Zをわずかに変えてもln(-ln(Z(t)))がさほど変化しないので、そのデータは軽視されることになる。もちろん、本当にこんな尺度で残差を測るのが目的にかなっているんなら、これで良い訳です。が、多分そうじゃないでしょう。

　さて、(a)(b)はどちらも非線形最小二乗法の問題であり、「E(もしくはD)を最小にするZのパラメータ」を近似値から出発して、E(もしくはD)がもっと小さくなるようにパラメータを少しずつ改良する繰り返し計算が必要です。幸いなことに(c)によってかなり正確な近似値が分かっていますから、少ない繰り返し回数で簡単に収束するでしょう。

　とは言っても「そんな大変そうな計算はやりたくない。線形回帰だけでなんとかしたい」というご注文でしたら、たとえば(a)の代わりに
(e)残差
　　γ(t) = ln(-ln(Y(t)))-ln(-ln(Z(t)))
の重み付き二乗和
　　H=Σ(w(t)γ(t))^2
を最小化する線形回帰を行う方法があります。ここで重みw(t)は
　　1/w(t) = (d/dY)(ln(-ln(Y)))　　(Y=Y(t))
と決めます。こうすれば、Y（あるいはZ）のわずかな変化がln(-ln( ))という曲線によって拡大・縮小される効果を概ね打ち消すことができます。なので、(c)よりも(a)にかなり近い結果が得られるでしょう。

お礼

リクエストにお答えいただきありがとうございます。

非線形最小二乗法の問題を繰り返し計算で解くことはハードルが高く、(e)の方法を検討したいと思っています。
計算をしたところ、w(t)=Yln(Y) であることはわかったのですが、どう線形回帰を使えばいいのかがわかりません。
「二乗和の微分=0 を解く」といったことを自力でやる必要がありますか？
可能なら、Excelの線形回帰の関数slope, intercept等で実現可能な方法が望ましいのですが。

よろしくお願いいたします。

回答No.1

stomachman

ワイブル分布の密度関数をフィッティングしてみたらなかなかうまいことあてはまり、m>1 になったと。ワイブル分布に従うならある時点tに破棄される確率密度は
m(η^(-m))(t^(m-1))
である。今後もワイブル分布に従うかどうかは、「破棄される確率が時間tのm-1乗に比例する」という定量的関係が成立つかどうかに掛かっています。
　さて、m>1ということは「この製品は、時間が経つにつれて破棄される確率が大きくなる傾向がある」ということを意味する。これは、たとえば製品が耐久消費財である場合には納得できる話です。どんな耐久消費財であれ、その製品が壊れるとか、もっと良い新製品がより安く手に入るようになるとか、単に飽きるとか、引っ越しをする際の大掃除だとか、持ち主が死ぬとか、破棄する理由の多くは時間が経つにつれて強くなっていくでしょうから。ですから、定性的には辻褄が合っています。
　定量的にはどうか。「ある年までいくつ残っているか」は、「前の年に残っていた個数」×「ある年に残っているもののうち何％が減るか（上記の式）」で計算されます。「毎年残っているもののうちの何％が減るか」の予測誤差は年々大きくなるでしょう。しかし、これと掛け算する「前の年に残っていた個数」が年々小さくなるために、誤差の影響が小さく抑えられる。つまり、「破棄される確率が時間tのm-1乗に比例する」ということが現実に対する近似としてあまり良い精度を持っていなかったとしても、残っている数がある程度減ってしまったあとなら、ソレナリに予測に使えそうですね。
　ただ、予測の目的にもよります。時間が経つにつれて破棄する理由が却って小さくなる要因としては「骨董的価値（ヴィンテージ）」というものがありますよね。製品によっては、これのおかげで、一定の年数が経過し残りがわずかになってくると、それ以上はなかなか減らなくなるでしょう。ですから、何十年何百年後にいくつ残っているか、ということまで予測しようと思ったら、ワイブル分布じゃ話にならないわけです。

　以上の議論から、お考えになったやり方は「今後も同程度の精度が出ると保証される訳ではないが、ま、使い途によってはソレナリにソレッポイ」という程度の意味があるでしょう。

　きちんとやりたければ、（過去の経験に基づく）統計のアプローチとして、もっと沢山の類似製品の長期に渡るデータを使って、経験則がどこまで成立ちそうかを見極めることが必要である。また、（過去も未来もひっくるめて考えられる）確率論のアプローチとして、現実をうまく説明するような「破棄される確率が時間tとともにどう変わるか」の理論モデルを考案して、ワイブル分布よりも精度の良い予測理論を検討する。そしたら立派な経済学の論文が書ける訳で、これは大変な研究になっちゃうでしょう。

お礼

お礼が遅くなり申し訳ありません。

ありがとうございます。

いただいたご指摘はごもっともで、「この予測方法は妥当だろうか」という問いに真正面からお答えいただき感謝してます。

質問の仕方が悪かったのですが、お聞きしたかったのは、「ワイブル分布をフィッティングする」ということは確定だとして、その具体的やり方が箇条書きで書いた方法であっているのか、ということでした。「回帰分析でいいの？」というのが一番のポイントです。
この視点でもお答えいただけたら幸いです。