2010-10-31

統計解析ツール ~ 歩き始め ~

元配置の分散分析ができたからといって、全て見通しが立ったわけではありません。しかし、自分が納得できる統計解析ツール、特に実験計画法の為の解析ツールを作りたいという欲求を抑えることができず、少しずつ GUI の部分を作り始めています。

GUI では、かゆいところに手が届くような操作性と柔軟性を実現したいので、じっくり時間をかけて作り込む予定です。若いころに BBN RS/1 というツールを VAX Station 上で使っていましたが、実験計画用のモジュール (Discover とか Explorer という名前でした)の使い心地が忘れられません。今なら GUI でもっともっと使い易く出来るはずです。



実は、むかしむかしも同じようなツールを自分で作りたいと考えていました。その時は応答曲面法 (RSM) を全面に出していましたが、行列計算の量が多く、手許の PC で実用的な処理速度を得られず、挫折してしまいました。

そもそも Tcl/Tk を使い続けている理由も、こういうツールを自分で作りたいという欲求があるからです。かつての人気はすっかり無くなってしまいましたが、Tcl/Tk の開発は止まっていません。Tktable も BLT (RBC) も利用できます。Eclipse などの統合開発環境が利用出来るので、ソフトウェアの開発が格段に楽になってきてもいます。

今回は、オーソドックスに一元配置の分散分析からはじめて、直交表までをカバーして、その次に中心複合計画やコンピュータの計算応力を駆使する最適計画など応答曲面法までを目指したいと考えています。

さて、完成はいつになるやら…。まずはアルファ版を公開することを目指します。とりあえずは決意表明をしておこう。

参考/引用文献

[1] 岩崎 学 著, 統計的データ解析入門 実験計画法 (東京図書)ISBN4-489-00725-6
 

2010-10-16

RBC - Refactored BLT Component

Fedora では Tcl/Tk 8.5 で BLT パッケージ(2.4z にパッチをあてたもの)が利用できますが、nightly-cvs の Tcl/Tk8.6 を自分でコンパイルして使用している環境では、同サイトで入手できる blt (3.0) をどうしてもうまくビルドできません。なんとかならないかとインターネット上で情報をあさっていたところ、古い記事ですが BLT 2.4z を リファクタリングした RBC についての議論を見つけました[1]。どうやら Mayo SPPDG - Mayo Clinic and FoundationRobert W Techentin 氏が、院生のチームと一緒になって取り組んだ成果のようです。

SourceForge.net にサイト[2] が開設されているので、早速ダウンロードしてビルドしてみました。試行錯誤の後、結局、以下のようにしてビルドしました。ビルドした環境は Fedora 14 (beta) x86_64 です。
$ tar zxvf rbc0.1.src.tar.gz
$ cd rbc
$ CFLAGS="-DUSE_INTERP_RESULT" \
> ./configure --prefix=$HOME --with-tcl=../tcl/unix --with-tk=../tk/unix
$ make
$ make install
コンパイルフラグ USE_INTERP_RESULT は、Tcl 8.6 ではレガシーな (Tcl_Interp) interp->result にアクセスする際に必要になります[3]
テスト用に BLT 用の以前用意したサンプルを RBC 用に書き換えて実行してみました。
#!/bin/sh
# the next line restarts using wish \
exec wish8.6 "$0" "$@"
 
package require rbc

proc init_x {} {
    set xList {}
    for {set i 0} {$i <= 200} {incr i} {
        lappend xList [expr $i / 5.]
    }
    return $xList
}

proc init_y {xList} {
    set yList {}
    foreach x $xList {
        lappend yList [expr sqrt($x) * sin($x)]
    }
    return $yList
}
 
# ----------------------------------------------------------------------------
# メイン
# ----------------------------------------------------------------------------
wm title . "rbc::graph"

set grp ".grp"
rbc::graph $grp \
        -title           "Sample Graph" \
        -font            "helvetica 16" \
        -rightmargin     15 \
        -plotbackground  "light cyan" \
        -plotpadx        0 \
        -plotpady        0 \
        -plotborderwidth 0
pack .grp

# ズームスタックを有効にする
#Blt_ZoomStack $grp

# X と Y 軸
$grp axis configure x \
        -subdivisions 10 \
        -tickfont     "courier 8" \
        -title        "X Axis" \
        -titlefont    "helvetica 12"
$grp axis configure y \
        -subdivisions 10 \
        -tickfont     "courier 8" \
        -title        "Y Axis" \
        -titlefont    "helvetica 12"

# グリッド
$grp grid configure \
        -color "turquoise" \
        -mapx  x \
        -mapy  y
$grp grid on

# 凡例
$grp legend configure \
        -font     "helvetica 10" \
        -ipadx    10 \
        -position bottom \
        -relief   ridge

# データ
set xList [init_x]
set yList [init_y $xList]
$grp element create func \
        -label        "sqrt(x) * sin(x)" \
        -xdata        $xList \
        -ydata        $yList \
        -smooth       quadratic \
        -linewidth    1 \
        -color        blue \
        -symbol       circle \
        -pixels       0.0625i \
        -scalesymbols yes \
        -outline      red \
        -fill         green

# ---
# sample_graph.tcl
実行例を以下に示しました。日本語フォントが表示できないし、 BLT のすべての機能が使えるわけでもなく制限はありますが、Tcl/Tk 8.6 の環境で手軽に(そこそこの)グラフを作成できそうです。

参考サイト

[1] RBC - Refactored BLT Components - Initial Release - Rhinocerus
[2] The RBC Toolkit | Download The RBC Toolkit software for free at SourceForge.net
[3] DIRECT ACCESS TO INTERP->RESULT - Tcl_SetResult manual page
 

2010-10-10

ちょっと分散分析(2)

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 $

# サンプルデータ1
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}]

# サンプルデータ2
#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 などの配布ライセンスで良いライブラリがあれば乗り換えようと思っています。オープンソースを推進したい者としては矛盾した考えですが、しかたがありません。
 

2010-10-02

ちょっと分散分析

々から作りたいと思っていて、なかなか重い腰を上げて取り組もうとしなかった分散分析のライブラリについて、ようやく少しずつ作り始めました。

本当は 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 $

# サンプルデータ1
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}]

# サンプルデータ2
#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