ワイブル分布のパラメータの最尤推定     Last modified: Dec 07, 2004

目的

ワイブル分布のパラメータを最尤推定する。

使用法

weibull.par(x, epsilon=1e-7)

引数

x         データベクトル
epsilon   収束判定値(省略時には 1e-7 が仮定される)
max.loop  収束計算の上限回数(省略時には 500 回が仮定される)

ソース

大きさ n の標本 x1, x2, ..., xn から,α および m の最尤推定量を求めるには,
1/α = n / Σ xim
m = n / (Σ(xim * ln(xi))-Σln(xi))
を,1/α と m に関して反復法により求める。

インストールは,以下の 1 行をコピーし,R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/weibull_par.R", encoding="euc-jp")

# ワイブル分布のパラメータの最尤推定
weibull.par <- function(x,                                   # データベクトル
                        epsilon=1e-7,                           # 収束判定値
                        max.loop=500)                           # 収束計算の上限回数
{
        x <- x[!is.na(x)]                                    # 欠損値を持つケースを除く
        n <- length(x)                                               # データの個数
        m <- a <- 1                                               # 初期値を設定する
        error <- TRUE                                                # 収束したかどうかのフラッグ
        for (i in 1:max.loop) {                                 # 収束するまで繰り返し
                ao <- a                                              # 改善前の値
                mo <- m                                              # 改善前の値
                temp1 <- x^mo
                temp2 <- log(x)
                a <- n/sum(temp1)                            # 改善後の値
                m <- n/(a*sum(temp1*temp2)-sum(temp2))               # 改善後の値
                if (abs(a-ao) < epsilon &&
                    abs(m-mo) < epsilon) {                   # 差が収束判定値以下ならば,
                        error <- FALSE                               # 収束した
                        break
                }
        }
        if (error) {                                            # 収束しなかった場合 error = TRUE のまま
                warning("収束しませんでした")
        }
        scale <- (1/a)^(1/m)                                 # 尺度パラメータ
        return(c("shape"=m, "scale"=scale))                     # m は形状パラメータ
}


使用例

x <- rweibull(90, 1.5, scale=3)	# 90 個のワイブル乱数を発生させる
weibull.par(x)

出力結果例

> x <- rweibull(90, 1.5, scale=3)
> weibull.par(x)	# 使用されるデータが毎回違うので,表示される解も毎回違う
   shape    scale 
1.537048 3.162941 

・ 解説ページ
注:このプログラムで計算される尺度パラメータは,
分布関数が F(x) = 1 - exp(- (x/b)^a) と表現されるときの b である(R の *weibull 関数もこの流儀)。
分布関数の定義によっては F(x) = 1 - exp(- (x^a)/b) とされることがある(上の解説ページの流儀)。
そのような場合には,得られる scale の値を a 乗すればよい(a = shape)。


・ 直前のページへ戻る  ・ E-mail to Shigenobu AOKI ( @si.gunma-u.ac.jp )

Made with Macintosh