目的

Box-Cox 変換の，最適のλをシンプレックス法によって求める

使用法

Box.Cox.transformation2(x, loop = 500, epsilon = 1e-15, Alpha = 2, Beta = 0.5, Gamma = 2)

引数

x	データベクトル
loop	収束計算の許容ループ数
epsilon	収束判定値
Alpha	シンプレックス法のα
Beta	シンプレックス法のβ
Gamma	シンプレックス法のγ

ソース

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

# Box-Cox 変換の，最適のλをシンプレックス法によって求める
Box.Cox.transformation2 <- function( x,                              # データ
                                        loop = 500,                     # 収束計算の許容ループ数
                                        epsilon = 1e-15,                # 収束判定値
                                        Alpha = 2,                      # シンプレックス法のα
                                        Beta = 0.5,                     # シンプレックス法のβ
                                        Gamma = 2)                      # シンプレックス法のγ
{
        Box.Cox.sub <- function(lambda)                                      # ボックス・コックス変換の計算値
        {
                if (lambda == 0) {
                        w <- Gm*log(x)
                }
                else {
                        w <- (x^lambda-1)/(lambda*Gm^(lambda-1))
                }
                sd(w)
        }

        x <- x[!is.na(x)]                                            # 有効データのみを対象とする
        Gm <- exp(mean(log(x)))                                              # 幾何平均を求める

        p1 <- -3                                                     # 初期値
        p2 <- p1+0.1                                                 # 初期値

        vec <- c(p1, p2)

        for (i in 1:loop) {                                             # 収束計算

                result <- c(Box.Cox.sub(vec[1]), Box.Cox.sub(vec[2]))

                h <- which.max(result)                                       # 大きい方
                s <- which.min(result)                                       # 小さい方

                ph <- vec[h]
                fh <- result[h]
                ps <- vec[s]
                fs <- result[s]

                p0 <- vec[s]

                pr <- (1+Alpha)*p0-Alpha*ph
                fr <- Box.Cox.sub(pr)

                if (fr > fh) {                                          # 探索値の状況により推測値を修正
                        pc <- Beta*ph+(1-Beta)*p0
                        vec[h] <- pc
                }
                else if (fs > fr) {
                        pe <- Gamma*pr+(1-Gamma)*p0
                        fe <- Box.Cox.sub(pe)
                        if (fr > fe) vec[h] <- pe else vec[h] <- pr
                }
                else {
                        vec[h] <- pr
                }
                if (abs((vec[1]-vec[2])/vec[1]) < epsilon) break     # 近似値の差が小さいなら収束と見なす
        }

        mean(vec)                                                       # 2つのベクトルの平均値は，真値により近いだろう
}


使用例

> x <- c(5.0, 5.0, 3.3, 4.3, 4.0, 5.5, 4.0, 6.0, 5.0, 5.0, 4.0, 4.3, 5.3, 5.0, 6.0, 6.7, 6.5, 6.0, 6.0, 5.3, 7.0, 5.0, 6.3, 5.3, 4.5, 6.0, 7.0,
+ 2.0, 2.5, 1.5, 1.7, 1.0, 1.0, 2.0, 1.0, 1.7, 2.0, 2.0, 1.0, 1.3, 2.5, 1.0, 2.0, 2.0, 1.0, 3.0, 1.3, 1.0, 2.0, 2.0, 1.7, 1.0, 1.0, 1.0,
+ 4.0, 5.3, 3.0, 3.7, 2.7, 3.3, 5.0, 2.7, 4.7, 3.7, 3.7, 4.0, 4.7, 4.3, 5.7, 4.7, 5.3, 5.5, 4.3, 7.0, 5.0, 5.3, 5.7, 5.3, 4.3, 5.0, 5.0, 5.0, 4.3, 4.7,
+ 1.3, 1.3, 1.3, 1.7, 2.0, 1.0, 1.3, 1.3, 1.7, 1.7, 1.3, 1.3, 1.3, 2.3, 2.0, 2.7, 2.0, 1.7, 2.3, 1.0, 2.0, 2.0, 2.0, 1.7, 2.0, 1.3, 1.3, 2.0, 1.3, 1.7,
+ 5.0, 5.0, 3.7, 3.3, 3.0, 3.5, 3.3, 3.0, 3.0, 4.3, 3.7, 4.0, 4.0, 3.3, 5.0, 6.0, 4.0, 4.0, 5.7, 5.5, 5.5, 5.5, 5.3, 4.7, 5.0, 4.3, 4.3, 3.0, 6.0,
+ 2.3, 1.7, 2.0, 1.0, 1.5, 1.0, 1.0, 2.0, 1.0, 2.0, 1.3, 1.7, 1.7, 2.0, 2.0, 2.0, 1.7, 2.0, 1.7, 0.5, 1.7, 1.3, 1.7, 2.3, 2.0, 1.0, 1.0, 1.3, 2.0,
+ 2.0, 5.5, 3.5, 5.3, 5.0, 6.0, 8.0, 6.0, 4.0,
+ 2.0, 3.0, 2.5, 2.0, 1.0, 2.0, 1.0, 3.0, 2.0)

> Box.Cox.transformation2(x)

[1] 0.2368628	# λの推定値

 注：得られたλをそのまま採用するのではなく，意味の明らかな近傍の値を採用するべきである。
　　　λ = 0 は対数変換
　　　λ = 0.5 のときは平方根変換
　　　λ = -1 のときは逆数変換となる
　　　λ = 1 のときは，単なる線形変換（もとのデータから 1 引くだけ）なので，分布型は変わらない。