前々から作りたいと思っていて、なかなか重い腰を上げて取り組もうとしなかった分散分析のライブラリについて、ようやく少しずつ作り始めました。本当は C# で作りたかったのですが、GUI 化の手間と、ゆくゆくは超越関数（Γ関数など）を扱う統計値の計算ルーチンを作ることを考えれば、慣れている Tcl/Tk (+ Tktable) と C の組み合わせで作る方が現実的と考えました。
また、TclOO を利用したかったので、Tcl の nightly-cvs のサイトから Tcl/Tk の最新のソース (8.6b1.2) をダウンロード、Fedora 上でビルドして利用しています。
まずは一元配置の簡単なルーチンを作ってみました。参考サイト [1] と [2] で検算をしています。
#!/bin/sh
# the next line restarts using tclsh \
exec tclsh8.6 "$0" "$@"
# 分散分析
# $Id: anova_1way.tcl,v 1.4 2010/10/02 15:00:00 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
# 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]
# 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
puts "ANOVA (1way)"
puts "Source SS df MS F 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 ""
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
}
# -------------------------------------------------------------------------
# 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 Between 0.318 2 0.159 4.615 Within 0.310 9 0.034 Total 0.628 11まだ、p値の計算はできていません。この部分は C で Tcl のライブラリを作成する予定です。あとは、すこしずつ GUI を作りながら、二元配置、直交表の分散分析にも対応させていく予定です。
参考サイト
[1] 一元配置分散分析[2] One Way ANOVA
0 件のコメント:
コメントを投稿