p 値の計算をするためにガンマ関数やらベータ関数を計算するルーチンを C で作成し、Tcl の拡張パッケージ上の拡張された算術関数として利用できるようにして、その上で P 値を計算するルーチンを作ってはみましたが、計算した値が合いません。きっと自分の理解不足によるもののと考え、とりあえず諦めてしまいました（泣）。
そのかわり近似計算をしているサイトをいろいろ探し、サイト[1] を参考にして P 値を与えるルーチンを追加しました。悔しいのでもっと勉強して、自分で納得できるルーチンを作成していきます。
#!/bin/sh
# the next line restarts using tclsh \
exec tclsh8.6 "$0" "$@"
# 分散分析
# $Id: anova_1way.tcl,v 1.6 2010/10/10 05:15:34 bitwalk Exp bitwalk $
# サンプルデータ１
set a [list {3.42 3.84 3.96 3.76}]
set b [list {3.17 3.63 3.47 3.44 3.39}]
set c [list {3.64 3.72 3.91}]
# サンプルデータ２
#set a [list {25 30 20 32}]
#set b [list {30 33 29 40 36}]
#set c [list {32 39 35 41 44}]
# データ
set data [list $a $b $c]
# -----------------------------------------------------------------------------
# クラス anova
# -----------------------------------------------------------------------------
oo::class create anova {
constructor {a} {
my variable x way
# データリスト
set x $a
# 分散分析の対象フラグのリセット
set way ""
}
# -------------------------------------------------------------------------
# 一元配置
# -------------------------------------------------------------------------
method oneway {} {
my variable x way
my variable n k
my variable SSb SSw SSt
my variable df_b df_w df_t
my variable MSb MSw
my variable F p
# Sum of Squared Deviations
# SSt (TOTAL)
set x_simple [my list_simple $x]
set n [llength $x_simple]
set x_mean [my mean $x_simple $n]
set SSt [my var $x_simple $x_mean]
# SSb (BETWEEN GROUPS)
set k [llength $x]
set SSb 0
for {set i 0} {$i < $k} {incr i} {
set x_b($i) [my list_simple [lindex $x $i]]
set n_b($i) [llength $x_b($i)]
set x_mean_b($i) [my mean $x_b($i) $n_b($i)]
set SSb [expr $SSb + $n_b($i) \
* ($x_mean_b($i) - $x_mean) \
* ($x_mean_b($i) - $x_mean)]
}
# SSw (WITHIN GROUPS)
# set SSw [expr $SSt - $SSb]
set SSw 0
for {set i 0} {$i < $k} {incr i} {
set x_w($i) [my list_simple [lindex $x $i]]
set n_w($i) [llength $x_w($i)]
set x_mean_w($i) [my mean $x_w($i) $n_w($i)]
set SSw [expr $SSw + [my var $x_w($i) $x_mean_w($i)]]
}
# Degree of Freedom
set df_b [expr $k - 1]
set df_w [expr $n - $k]
set df_t [expr $n - 1]
# Mean Square
set MSb [expr $SSb / $df_b]
set MSw [expr $SSw / $df_w]
# F value, the ratio of the two variance estimates
set F [expr $MSb / $MSw]
set p [my f_pval $F $df_b $df_w]
# flag
set way "One Way"
return 1
}
# -------------------------------------------------------------------------
# result
# 結果の表示
# -------------------------------------------------------------------------
method result {} {
my variable way
switch $way {
"One Way" {
my result_1way
}
default {
puts "No result yet!"
}
}
}
# -------------------------------------------------------------------------
# result_1way
# 結果の表示（一元配置）
# -------------------------------------------------------------------------
method result_1way {} {
my variable x
my variable n k
my variable SSb SSw SSt
my variable df_b df_w df_t
my variable MSb MSw
my variable F p
puts "ANOVA (1way)"
puts "Source SS df MS F val P Val"
puts -nonewline "Between"
puts -nonewline [format " %8.3f" $SSb]
puts -nonewline [format " %8d" $df_b]
puts -nonewline [format " %8.3f" $MSb]
puts -nonewline [format " %8.3f" $F]
puts -nonewline [format " %8.6f" $p]
if {$p < 0.001} {
puts " ***"
} elseif {$p < 0.01} {
puts " **"
} elseif {$p < 0.05} {
puts " *"
} else {
puts ""
}
puts -nonewline "Within "
puts -nonewline [format " %8.3f" $SSw]
puts -nonewline [format " %8d" $df_w]
puts -nonewline [format " %8.3f" $MSw]
puts ""
puts -nonewline "Total "
puts -nonewline [format " %8.3f" $SSt]
puts -nonewline [format " %8d" $df_t]
puts ""
}
# -------------------------------------------------------------------------
# mean
# 平均値の計算
# -------------------------------------------------------------------------
method mean {xList num} {
set sum [my sum $xList]
return [expr $sum / $num]
}
# -------------------------------------------------------------------------
# sum
# 合計の計算
# -------------------------------------------------------------------------
method sum {xList} {
set total 0.0
foreach x $xList {
set total [expr $total + $x]
}
return $total
}
# -------------------------------------------------------------------------
# var
# 分散の計算
# -------------------------------------------------------------------------
method var {xList xMean} {
set var 0.0
foreach x $xList {
set var [expr $var + ($x - $xMean) * ($x - $xMean)]
}
return $var
}
# -------------------------------------------------------------------------
# f_pval
# P値（F分布の上側確率の算出）
# -------------------------------------------------------------------------
method f_pval {f df1 df2} {
set m2pi 0.636619772367581343076
set m_pi 3.14159265358979323846
set i1 [expr 2 - ($df1%2)]
set i2 [expr 2 - ($df2%2)]
set w [expr $df2/($df2 + $df1*$f)]
set y [expr 1.0 - $w]
set p [expr sqrt($w)]
set s [expr sqrt($y)]
if {[expr 2*$i1 - $i2] == 0} {
set p [expr 1.0 - $s]
set c [$s*$w/2.0]
} elseif {[expr 2*$i1 - $i2] == 1} {
set c [expr $s*$p/$m_pi]
set p [expr 1.0 - atan2($s, $p) * $m2pi]
} elseif {[expr 2*$i1 - $i2] == 2} {
set p $w
set c [expr $w*$y]
} else {
set c [expr $y*$p/2.0]
}
for {set i2 $i2} {$df2 > $i2} {incr i2 2} {
set p [expr $p - 2.0/$i2*$c]
set c [expr $c*$w*($i1 + $i2)/$i2]
}
for {set i1 $i1} {$df1 > $i1} {incr i1 2} {
set p [expr $p + 2.0/$i1*$c]
set c [expr $c*$y*($i1 + $i2)/$i1]
}
return [expr ($p < 0.0 && abs($p) < 1e-10) ? 0.0 : $p]
}
# -------------------------------------------------------------------------
# list_simple
# 入れ子になっているリストをシンプルなリストに変換
# -------------------------------------------------------------------------
method list_simple {xList} {
regsub -all {[\{\}]} $xList {} foo
regsub -all {\s{2,}} $foo { } baa
return $baa
}
}
# -----------------------------------------------------------------------------
# メイン
# -----------------------------------------------------------------------------
set obj [anova new $data]
if {[$obj oneway]} {
$obj result
} else {
puts "failed"
}
# ---
# END PROGRAM
実行例を以下に示しました。
$ ./anova_1way.tcl
ANOVA (1way)
Source SS df MS F val P Val
Between 0.318 2 0.159 4.615 0.041749 *
Within 0.310 9 0.034
Total 0.628 11
参考サイト
[1]
統計関数の確率計算
[2]
The Free Statistics Calculators Website - Home （算出値の検証用）
追記
この記事を書いた後、なぜ用意した関数でうまくできなかったのかが突然判ってしまいました。長年 Tcl を使っているのに、Tcl 特有の計算の落とし穴にすっかりハマっていたのです。
以下は、用意していた規格化された不完全ベータ関数
betai を使って書き直したメソッド
f_pval です。Tcl では、整数÷整数の結果は自動的に整数になってしまいます。浮動小数の結果が必要な場合は、分母分子のどちらか一方を浮動小数の値に明示にしておく必要があります。
# -------------------------------------------------------------------------
# f_pval
# P値（F分布の上側確率の算出）
# -------------------------------------------------------------------------
method f_pval {f df1 df2} {
set x [expr $df1 * $f / ($df2 + $df1*$f)]
return [expr 1 - betai($df1/2.0, $df2/2.0, $x)]
}
実行例を以下に示しましたが、上記と同じ結果が得られています。
$ ./anova_1way.tcl
ANOVA (1way)
Source SS df MS F val P Val
Between 0.318 2 0.159 4.615 0.041749 *
Within 0.310 9 0.034
Total 0.628 11
今のところ規格化された不完全ベータ関数
betai の算出には、下記の
GSL の関数を利用しています。
double gsl_sf_beta_inc (double a, double b, double x)
以下のようなルーチンを使って Tcl の拡張パッケージを作成しています。
/*****************************************************************************
* int func_betai
* 規格化された不完全ベータ関数の処理部
*
* Tcl 側のコマンド
* expr betai($a, $b, $x)
*****************************************************************************/
int
func_betai (ClientData clientData, Tcl_Interp * interp,
Tcl_Value * args, Tcl_Value * resultPtr)
{
double a, b, x;
a = args->doubleValue;
args++;
b = args->doubleValue;
args++;
x = args->doubleValue;
/* 計算結果 */
resultPtr->type = TCL_DOUBLE;
resultPtr->doubleValue = gsl_sf_beta_inc (a, b, x);
return TCL_OK;
}
ただ、GSL の配布ライセンスは GPLv3 (
GNU General Public License) に従っているので、扱い難いと感じています。計算方法の検証ができたので、他に BSD などの配布ライセンスで良いライブラリがあれば乗り換えようと思っています。オープンソースを推進したい者としては矛盾した考えですが、しかたがありません。
にほんブログ村
