Box-Cox変換(χ2値による補正を行うバージョン)は、データ分布の正規性を増し、かつブロック間の等分散性も増加させたい場合に使います。 データ分布の正規性を増すことだけが目的でしたら(あまり、そういう場面は無いと思いますが)、これではなく通常のBox-Cox変換をお使い下さい。
# Box-Cox transformation with an adjustment of chi-square values # (C) Copyright 2002, Hisashi SATO # # B_test()は、Bartlett's χ二乗値を求めるメソッド def Box_Cox2(array, add) # 定数、変数の初期設定 length = 1.0 # ポイント間の初期距離 ramda_dn = -10.0 # λ値(下位置) ramda_md = ramda_dn + length # λ値(中位置) ramda_up = ramda_md + length # λ値(上位置) # 配列の準備と全サンプル数の数え上げ、自由度の計算 sample_number = 0 # 全サンプル数 transformed_dn = Array.new() # ramda_upで変換されたデータ値 transformed_md = Array.new() # ramda_mdで変換されたデータ値 transformed_up = Array.new() # ramda_dnで変換されたデータ値 for block in 0..array.size-1 transformed_dn << Array.new(array[block].size,0.0) transformed_md << Array.new(array[block].size,0.0) transformed_up << Array.new(array[block].size,0.0) sample_number += array[block].size end df = sample_number-array.size # 自由度 # データに0が含まれていたら、小さな数を足す for block in 0..array.size-1 for n in 0..array[block].size-1 if array[block][n] == 0 then array[block][n] += add end end end # ■↓■↓■↓■3点法ループの開始■↓■↓■↓■ until length < 0.0001 || ramda_md >= 10 # 変換値とその合計を算出 sum_dn, sum_md, sum_up, sum_ln = 0.0, 0.0, 0.0, 0.0 for block in 0..array.size-1 for n in 0..array[block].size-1 if ramda_dn != 0 then transformed_dn[block][n] = (array[block][n] **ramda_dn - 1) / ramda_dn else transformed_dn[block][n] = log(array[block][n]) end if ramda_md != 0 then transformed_md[block][n] = (array[block][n] **ramda_md - 1) / ramda_md else transformed_md[block][n] = log(array[block][n]) end if ramda_up != 0 then transformed_up[block][n] = (array[block][n] **ramda_up - 1) / ramda_up else transformed_up[block][n] = log(array[block][n]) end sum_dn += transformed_dn[block][n] sum_md += transformed_md[block][n] sum_up += transformed_up[block][n] sum_ln += log(array[block][n]) end end # 変換値の分散を算出 mean_dn = sum_dn / sample_number mean_md = sum_md / sample_number mean_up = sum_up / sample_number varianth_dn, varianth_md, varianth_up = 0.0, 0.0, 0.0 for block in 0..array.size-1 for n in 0..array[block].size-1 varianth_dn += (transformed_dn[block][n] - mean_dn)**2 varianth_md += (transformed_md[block][n] - mean_md)**2 varianth_up += (transformed_up[block][n] - mean_up)**2 end end varianth_up = varianth_up / (sample_number-1) varianth_md = varianth_md / (sample_number-1) varianth_dn = varianth_dn / (sample_number-1) # L値を算出(この部分だけ、通常のBox_Cox変換プログラムと異なる) l_dn = -0.5*B_test(transformed_dn) - 0.5*df*log(varianth_dn) + (ramda_dn-1)*df*sum_ln / sample_number l_md = -0.5*B_test(transformed_md) -0.5*df*log(varianth_md) + (ramda_md-1)*df*sum_ln / sample_number l_up = -0.5*B_test(transformed_up) -0.5*df*log(varianth_up) + (ramda_up-1)*df*sum_ln / sample_number # L値の評価と次のλ値の設定 if l_dn < l_md and l_md < l_up then ramda_dn += length ramda_md += length ramda_up += length elsif l_dn > l_md then length = length / 4 ramda_md = ramda_dn + length ramda_up = ramda_md + length else length = length / 2 ramda_md = ramda_dn + length ramda_up = ramda_md + length end end # ■↑■↑■↑■ループの終わり■↑■↑■↑■ return ramda_md end # 心得のある方は、ちょっと読んで頂ければ分かるのですが # このプログラムでは、一度計算したL値を再び計算し直したりしているので # 計算量が必要最低量の倍以上も多くなっています。 # λ値−L値のhashを使って、改善しようとも考えたのですが # 3点法ルーチン部のコードが、かなり読みにくくなることが予想されるので # 結局、そのままとなっています。