目的

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

使用法

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）。