delphiお よび excel と統計処理

各種分布のｐ値

統計処理をするとき、計算機でやるからには p値は計算機で求めたい。

岡本 安晴「Delphiで学ぶデータ分析法」ＣＱ出版社、1998
がまず参考になります。クリックして目次を見てください。
標準正規分布、カイ2乗分布、t分布、F分布、Uの分布（Mann-Whitney U test）、Tの分布（Wilcoxon signed ranks test）、2項分布の確率を求めることができます。www でも公開されました。
　以前，岡本さんの精度(1e-17)では収束しないということをここに書いていましたが，コプロの精度がなぜか低くなっていたためでした。これに対応するには
(1)精度を低くする。 1e-15 なら収束する
(2)procedure TForm1.OKButtonClick　のbegin の後に
Set8087CW($1332);
を入れることによって解決します。

delphi の計算精度を簡易チェックするプログラムをつくりました(delphi 6)．seidomain2.exe
button1 を押してください． 数値プロセッサ指定前の値と指定後の値を表示します．
私の環境のwindows98 で走らせると

original data	1.12345678901234567890
$1272	1.1234567890123456	指定前
$1372	1.12345678901234568	指定後

という結果になります．

正規分布（英語、日本語）　t分布解説（英語 日本語）
F分布解説（英語、日本語） χ2分布（英語、日本語）
英語の解説はThe CRC Concise Encyclopedia of Mathematicsのサイト(今は閉じられています)
日本語の解説は青木繁伸＠群馬大学のサイトにリンクしています。

ただし、この本の付録についている分布関数の収束規準が1e-17 になっていてこれでは収束しない（いつも忍耐力の限界を越えるので収束しないと言い切ることはできません）ので1e-15 に下げるとすぐに結果がでます。

岡本安晴氏の本付属のプログラムでも求めることができますが、コピーライトの主張が本にはありますから、そのまま公開可能かというと普通はダメです。そこでフリーのソフトが欲しい。

そこで次のサイトのソフトです。

TP Math (Mathematical library for Turbo Pascal)

この数値計算ライブラリは充実している。しかもfreeware宣言している。

統計の分布計算も多く含んでいる。ただしturbo pascal ですが、そのままで分布計算をすることも可能です。正規分布,ｔ分布、Ｆ分布、カイ２乗分布などがあるので代表的なものはこれですむ。

分布計算は mathspec.inc にある。これを使うには、mathcons.inc の定数と math387.inc の最後にある

function MathError : Integer;

procedure InitGlobalVar;

を追加してやる。intpowという関数があるので unit math を加えている場合はこれを power に書き換えてやる。もしmath がないなら(standard version) 数学関数ルーティン(math387.inc)も一緒に組み込んだ方がいいでしょう。

最初にInitGlobalVar を走らせておく。

このルーティンを使ってp値を求める delphi アプリをつくってみました。
ダウンロード（pvalue.lzh）(1998/10/18)　自分でコンパイルするのには別途mathspec.inc が必要です。実行ファイルはついていますので、使うだけならコンパイルの必要はありません。
(1998/10/19)distribution.pas を付け加えたので、自分でコンパイルする場合は、mathspec.inc を同一フォルダーに入れるだけでいいようになりました。
正規分布両側確率、ｔ分布両側確率、Ｆ分布上側確率、カイ２乗分布上側確率。mathspec.inc にあるポアソン分布、２項分布も求めることは可能。

mathspec.inc に入っている関数一覧

{ ************************** Special functions ************************* }

function Fact(N : Integer) : Float; { Factorial }

function Binomial(N, K : Integer) : Float; { Binomial coef. C(N,K) }

function Gamma(X : Float) : Float; { Gamma function }

function SgnGamma(X : Float) : Integer; { Sign of Gamma function }

function LnGamma(X : Float) : Float; { Log(|Gamma(X)|) }

function IGamma(A, X : Float) : Float; { Incomplete Gamma function }

function JGamma(A, X : Float) : Float; { Complement of IGamma }

function Beta(X, Y : Float) : Float; { Beta function }

function IBeta(A, B, X : Float) : Float; { Incomplete Beta function }

function Erf(X : Float) : Float; { Error function }

function Erfc(X : Float) : Float; { Complement of Erf }

{ ********************** Probability distributions ********************* }

{ Binomial distribution with probability P and number of repetitions N }

function PBinom(N : Integer; P : Float; K : Integer) : Float; { Prob(X = K) }

function FBinom(N : Integer; P : Float; K : Integer) : Float; { Prob(X ＜= K) }

{ Poisson distribution with mean Mu }

function PPoisson(Mu : Float; K : Integer) : Float; { Prob(X = K) }

function FPoisson(Mu : Float; K : Integer) : Float; { Prob(X ＜= K) }

{ Standard normal distribution }

function FNorm(X : Float) : Float; { Prob(U ＜= X) }

function PNorm(X : Float) : Float; { Prob(|U| >= |X|) }

{ Student distribution with Nu d.o.f. }

function PStudent(Nu : Integer; X : Float) : Float; { Prob(|t| >= |X|) }

{ Khi-2 distribution with Nu d.o.f. }

function PKhi2(Nu : Integer; X : Float) : Float; { Prob(Khi2 >= X) }

{ F distribution with Nu1 and Nu2 d.o.f. }

function PSnedecor(Nu1, Nu2 : Integer; X : Float) : Float; { Prob(F >= X) }

カイ二乗分布上側確率を求める別のアルゴリズムのプログラムも作っています。青木繁伸さん（群馬大学のawk プログラムをdelphi に移植しました。ダウンロード(1998/8/25)
delphiにはこのタイプでの確率計算をするプログラムがなかなか見つかりません。サブルーティンとして組み込むためのものですが、一応カイ２乗値と自由度を入れると確率(p値）を計算する実行プログラムをつけています。

統計分布関数の精度をチェックするには，
DSTATTAB 含む非心(dos,win32,mac)

Computation of Statistical Distributions (ELV) dos

Probability Calculator　NCSS 統計分布関数・逆関数計算機　非心，studentized range を含む　Windows 一番手軽

tetrachoric correlation

John Uebersax氏の The Tetrachoric and Polychoric Correlation Coefficients にある,TCORR.FOR というfortan プログラムを excel vba に変換しました。tetra_corr.bas

正規分布を前提とした相関係数です。説明はUebersax氏のサイトを読んでください。

Dave Garson 氏のCorrelation 説明も読むといいかな。こっちは簡単な説明。

tetrachoric correlation(四分相関)行列の因子分析については，Factor Analysis of Tetrachoric Correlations

pascal のプログラムも公開されている。Dirk Enzmann　のサイトの一番下。文献については、SPSSのマクロのr_tetra.sps の中に複数あがっている。頻度が５未満のセルがあるとき精度が落ちるということだ。このpascal プログラムは対策済みのもの。実行ファイルは日本語モードでは動かないようだ。ソースファイルがあるのでそれをdelphi 化できる。delphi化した実行ファイルのみを参考出展。処理できる変数数を100変数まで拡張。データは1とそれ以外の２分をする。0,1のコーディングでも1,2のコーディングにも対応する。ただし、欠損値のあるデータは除いておくこと。 使用法は起動後、データファイルを指定。データファイルはスペースで区切っているデータのみのファイルであること。

polychoric correlation　多分相関係数(順序のあるクロス表の相関係数)

上のと同じくJohn Uebersax氏の The Tetrachoric and Polychoric Correlation Coefficients にある,XPC.FOR というfortan プログラムを excel vba に変換しました。polycorr1.xls(excel97 file) 個別に落とすならこっち（polycorr11.bas polycorr1.frm polycorr1.frx optionメニュー が必要ないならpolycorr11.bas だけでいいです。）

正規分布を前提とした順序尺度の相関係数です。２段階評定（２件法），３段階評定（３件法）や４段階評定（４件法）の尺度間の相関を求めます。説明はUebersax氏のサイトを読んでください。

Dave Garson 氏のCorrelation 説明も読むといいかな。こっちは簡単な説明。

おっと，狩野裕さん＠大阪大学の[fpr 294] FA of bin data に利用についていい説明があった。

tetrachoric correlation(四分相関)行列の因子分析についての説明，Factor Analysis of Tetrachoric Correlations 参照してください。
vba に移植してからいくつかチェックをいれた。
(1)２つの推定法のML（最尤）法と最小χ²法ではχ²法の推定値に納得いかない推定値があり，ML法を使うのがいいであろう。
(2)2-step 推定法とjoint 推定法ではjoint 推定法は収束がはっきりと遅い場合が多くある。１０倍以上時間がかかったりする。2-step推定法がいいであろう。
(3)cdf のルーティンとして，Uebersax氏が最初から入れているのは３つであるが，１つは落とした。alnorm, NPROBをそのまま入れている。ほかにexcel のnormsdist を追加した。excel のnormsdist は最速というわけでもない。alnormは精度に少し問題がある。NPROBを使うのがいいであろう。
(4)χ²値のp値を求めるルーティンはexcel のchidist に差し替えている。
(5)メタセル，初期値等は使用できない。書き出しについても必要があればもとのfortran プログラムを参考にして追加してください。
(6)このプログラムをpascal 化するのは骨でしょう。basic はfortran と親和性があるので比較的簡単であったが苦労はした。

使用法：
(a)相関係数だけを求めるときは関数 polycorr を使うと求めることができる。=polycorr(範囲指定)
(b)いろいろな指標も求める場合はマクロ, poly_corr を起動して範囲を指定する。
(c)optionを一時的に変えたいときはマクロ,poly_corr_option を起動してオプションを指定し，範囲を指定する。
(d)option を変えてしまいたいときは sub init0 のなかの変数の値を書き換える。
範囲はいずれもクロス表を指定する。マウスを使うこともできる。

poly_corr_Matrix_from_raw_data （polycorr11.bas のみに含まれている）

素データ行列から多分相関係数行列を求める。
(a)データ以外のセルからマクロ poly_corr_Matrix_from_raw_data を起動する。
(b)データは変数名が最初についていればそれを使うなければ，勝手に変数名をつける。
(c)データ範囲の指定はデータの下および右がempty ならば左上のセルを指定すればいい。
(d)データの下がempty なら，一列（変数名またはデータ）を指定すればいい。
(e)すべてのデータの範囲を指定してもいい。
(f)データに欠損値がないようにあらかじめ処理しておくこと。
(クロス表作成部分は青木繁伸氏＠群馬大学のvba プログラムを一部改変使用しています。)

prelis の公表している結果と比較すると，上記の方法と結果が一致した。ただし，cdf のルーティンはどれでも対して変わらない。

EQS の polychoric correlationの計算結果とは一致しないが，かなり近い値になる。

モンテカルロ法でpolychoric も求めるｃ言語プログラムが公表されている。

exact test

正確テストとか厳密検定とか直接確率計算法とか言われている。

Fisher の正確確率検定（直接法）(1997/9)　説明

青木繁伸さん（群馬大学）のAwk プログラムを移植しました。
一般に正確確率検定は片側検定になっていることが多いのですが、青木氏のプログラムは両側検定です。また、２×２だけでなく、ｍ×ｎに拡張されています。

多重比較

Tukey のスチューデント化された範囲の検定をfortran から移植しました。
ダウンロード（tukey.lzh）(1998/10/18 exeファイルを入れ忘れていたので追加しました。 2000/1/12)
多重比較の場合にもっともよく利用されるのが Tukey のスチューデント化された範囲 q　です。
Copenhaver,M.D. and Holland,B. 1988 Computation of the distribution of the maximum studentized range statistic with application to multiple significance testing of simple effects. J. Statist. Comput. Simul., 30, 1-15.
のfortranプログラムがStatLib（http://lib.stat.cmu.edu/） general の qprob(http://lib.stat.cmu.edu/general/qprob)にある。また、StatLibの-Software and extensions for the S(Splus) language(http://lib.stat.cmu.edu/S/)のtukey(http://lib.stat.cmu.edu/S/tukey)はこれを使用したプログラムである。計算は４倍精度から倍精度に移行してあるが、データ文の数値はこちらには載っているのでこれを使用した。これをdelphi(pascal)に移植したのがqprob1.pas である。

fortran のqerf, qerfc については大浦 拓哉（京都大学数理解析研究所）氏のOoura's Mathematical Software Packagesの「ガンマ関数 / 誤差関数 / その他初等関数」(http://momonga.t.u-tokyo.ac.jp/~ooura/gamerf-j.html） fortranプログラムをpascal に移植した。これもqprob1.pas に含めている。

岡本 安晴氏もdelphi プログラムを公表している（下の方）。ここには、サンプルサイズが不揃いな場合のq値を求めるプログラムもある。一般には手にはいらないから貴重(下に高速バージョンがある)。

永田・吉田『統計的多重比較法の基礎』サイエンティスト社(1997)の式をそのままプログラムにしたものである。同書によるとサンプルサイズが著しく変わらない場合は通常の求め方でもたいして変わらない。

同書は日本では初めての多重比較の専門書であり、多重比較をする場合は読む必要があるでしょう。how to ものとしても使えるように配慮している。

多重比較をする場合の必要サンプル数を求めるfortranプログラムおよびwww統計がJason C. Hsu's Home Pageにある。
Hsu,J.C.(1996) Multiple Comparisons: Theory and methods.Chapman & Hall
のプログラムである。そのほか同書で使用したデータおよびfortranプログラムが公開されている(上の書名をクリック)。

ここで公開されているのは、(a)全ての一対比較(Tukey)、(b)最大効果群との比較(Hsu)、(c)コントロール群との比較(Dunnett法)の３つの方法についてのものである。

一元の場合だけでなく、２元の場合の多重比較についても同書は述べている。が結局、SASの使用法の説明になる。repeated についても少し記述がある。

２元のサンプル数不揃いの場合のDunnett の 多重比較法のパーセンテージ・ポイントを求めるfortranプログラム(tmcc)および２元のサンプル数不揃いの場合の全ての対比較をする多重比較法のパーセンテージ・ポイントを求めるfortranプログラム（tpmc)がStatlib/General Archiveに公開されている。

Dunnett(1989)のfortranプログラムは、Statlib/Applied Statistics/ に公開されている(251)。
Multivariate normal probability integrals with product correlation structure.
使い方の説明については、
Bechhofer,R.E., Santner,T.J. and Goldsman,D.M.(1995) Design and analysis of experiments for statistical selection, screening, and multimple comparisons. Wiley. p277-278.を参照。

Toothaker,L.E.(1992)Multiple comparison procedures. Sage.
前提条件の満たされない（サンプル数が不揃い、母分散が等しくない、正規分布でない）ときと頑健性についてきちんと記述している。また、２元分散分析、ランダムブロック、繰り返し測定の場合についてもわかりやすく説明している。

もちろん、Hochberg,Y. and Tamhane,A.C.(1987) Multiple compariosn procedure. Wiley.
はもっとも詳しい。詳しすぎるので読みにくいのが欠点。

繰り返し測定の場合については、
Keselman,H.J. and Keselman,J.C.(1993) Analysis of repeated mesurements. in Edwards, L. K.(1993) Applied Analysis of Variance in Behavioral Science.Dekker
の4.2.4 Multiple comparison procedure (p.111-114),4.3.3 Multiple comparison procedure (p119-122)が詳しい。

インターネットでもtutorial が公表された。A Bluffer's Guide to ... Sphericity の最後のほう。

Tukey の有名なメモ The problem of multiple comparison(1953) は
Braun,H.I.(ed.) (1994)The collected works of John W. Tukey. vol.8.: Multiple comparisons:1948-1983,Chapman & Hall.
に含まれている。この巻は多重比較関係の論文を集めている。中心は上記メモでこれだけで300ページもある。総ページ数445ページ。

非心分布

検定力を求めるときや必要サンプル数を求めるときには非心分布が必要となる。ただ、ここに載せている、そのパラメータのときの確率を求めるだけでなく、その確率のときのパラメータを求めることが必要である。そのときは、例えば、ある確率のstudent化されたｑ値を求めるプログラムが参考になる。

非心t分布(1998/10/18)

fortran プログラム ALGORITHM AS 243 APPL. STATIST. (1989), VOL.38, NO. 1からの移植

非心χ2分布のpascal プログラムはhttp://lib.stat.cmu.edu/apstat/の204,231 にある。delphi で使用するのには修正が必要。
231のプログラムの関数centnorm は204 にある procedure を書き換えるといい。正規分布の±z の内側確率を求める。

非心F分布はFortranのプログラムしか見ていない。非心F分布の式（下の方）
DCDFLIB：Section of Computer Science, Department of Biomathematics, University of Texas M.D. Anderson Hospital.
のすざましいfortran program を移植しなければならないのだろうか？　すぐ下にあるDSTATTABをチェックしたところ、cumfnc.f を移植するだけですんだ。TP Math (Mathematical library for Turbo Pascal) のmathspec.inc を使用して、function をmathspec.inc に合わせる(ibeta と LnGamma)。

　非心χ2分布もこれについている、cumchn.f を書き直せばいい。

cumfnc.f を excel vba に移植した。非心ｆ分布累積確率関数。=ncdff(F値, df1, df2, 非心度)
cumchn.f を excel vba に移植した。非心χ²分布累積確率関数。=ncdchi(chi2値, df, 非心度)
cumtnc.f を excel vba に移植した。非心ｔ分布累積確率関数 =ncdft(t値 , df , 非心度)

DSTATTABにはdos版とＭａｃ版の実行ソフトもある。

任意の相関行列の乱数を作成する(choleskey 分解と乱数生成) excel vba

vba program(Choleskey_decomposition.bas)
任意の相関行列の乱数データを作成するにはチョロスキー分解をしてそれを使って乱数データを合成します。コレスキー分解をするvba が見つからなかったのでFortranプログラムをvbaに移植しました。乱数生成ルーチンは 伏見正則(『乱数』東京大学出版会, 1989)のプログラムを使用した。
　相関行列をそのまま再現する乱数行列にするにはあらかじめ最初の乱数行列を主成分分析をかけて無相関の乱数行列にしておく必要がある。spssマクロはran.sps

parallel analysis 因子分析の因子数決定法

ダウンロード(vba 部分はお見せできません）

因子分析における因子数決定法としてMAP, parallel analysis が最良の２つというデータがある。そのうち、parallel analysis(以下 PA)をエクセルのvbaで組んだ。しかし、理由があって中はお見せしない。

処理は対角１を入れたPAと対角SMCを入れたPAの２つである。対角１を入れたPAは、とりあえずのもっとも因子数をさしていると考える。対角SMCのPAは、最大可能の因子数をさしていると考える。parallel analysis において書いたように、唯一の絶対的な推定法はないので、いくつかの推定法とスクリープロットの状態などによって因子数を推定することになる。

このプログラムでは、各種重回帰式によるPAの推定値も出力する。
対角１の場合、Allen & Hubbard(1986)、Lautenshlager et al(1989)、Longman et al.(1989)PA mean、Longman et al.(1989)95%点、Keeling(2000)。結構予測がバラバラで笑えます。Allen & Hubbard(1986)の予測精度をよくしたものがLautenshlager et al(1989)ですからAllen & Hubbard(1986)はすでに必要のないものです。

対角SMCの重回帰による予測はMontanelli et al(1976)を出力します。

テストしたいデータの変数の数とサンプル数およびデータ行列の作成数を入れて、対角１、対角SMCを選択してから実行ボタンを押す。既定値としてはデータ行列の作成数を５０、対角１を選択している。処理結果は{P.A.2}のようにP.A. のあとに数字が追加されたシートに出力される。

結果は、PAの平均値、95％点を出力する。95%点は平均値＋1.645*SD の式によって推測している。平均値は比較的少ないデータ行列作成数で安定するようである。処理時間は変数の数が増えると増大する。

正規分布の乱数の発生法は伏見(1989)に紹介されている、Marsagliaの極座標法である。Box-Muller に比べ多少速いとある。

固有値のリストを残しますので、参考にしてください。

出力されたもののうち９５％点のほうを因子数を決定したいデータの相関行列の対角１の固有値、対角SMCの固有値と比較し、PAの固有値のほうが大きくなる前の固有値までが因子として有効です。

対角SMCの固有値を出せないソフトもあります。例えばSPSS。その場合、マクロで作るなりしてください。

SPSS、SAS、MATLABのPAを求めるマクロをO'Conner 氏が作ってます。そちらのほうが計算速度が速い。氏の場合、PAの95%点はSDからでなく、モンテカルロした結果の分布から求めています。また、SPSSにおいては SMCを入れた固有値も求めることもできます。また、ＭＡＰのプログラムもありますので、これらのソフトを使っている人は是非氏のサイトにあるマクロを利用してみてください。

[文献]
伏見正則(1989).　乱数　東京大学出版会
その他の文献は http://www.ec.kagawa-u.ac.jp/~hori/yomimono/pa.html を参照してください。

多変量解析・多次元尺度解析

因子分析、数量化３類、４類のプログラムは岡本 安晴氏のホームページ

delphi & pascal link

岡本 安晴 のホームページ

因子分析、数量化３類、４類、ステューデント化された範囲 q、

いろいろなシミュレーション実験などいろいろ delphi プログラムがある。

TP Math (Mathematical library for Turbo Pascal)

この数値計算ライブラリは充実している。しかもfreeware宣言している。

統計の分布計算も多く含んでいる。

Generic Turbo Pascal programming language material turbo Pascal の数値計算など揃っている

excel link

Accuracy of Statistical Packages Knusel　excel97 統計分布関数の評価

この論文で指摘されていることは，私のbinom.bas や青木さんの正規分布， F分布の関数においてはクリアされている。excel 2000でクリアされているかどうかは知らない。

Microsoft Excel を使った統計解析 青木繁伸＠群馬大学

VBA によるユーザ関数とアプリケーションプログラムの例　 青木繁伸＠群馬大学

揃ってるよね
関数
1.正規分布，カイ二乗分布，t 分布，F 分布の確率，および，パーセント点
2.2 群の平均値の差の検定を行う関数
3.独立性の検定（カイ二乗検定）を行う関数
4.フィッシャーの正確確率検定（2×2分割表）を行う関数
5.スピアマンの順位相関係数を求める関数
6.ケンドールの順位相関係数を求める関数
7.複数の相関係数の同等性の検定と母相関係数の点推定を行う関数
8.固有値・固有ベクトルを求める関数　（プロシージャ版）
プログラム
1.2 群の代表値の差の検定（t 検定，U 検定）
2.一元配置分散分析
3.クロス集計と独立性の検定
4.重回帰分析（ステップワイズ変数選択）
5.多項式回帰
6.判別分析（ステップワイズ変数選択，線形判別関数）
7.主成分分析
8.因子分析
9.数量化 I 類
10.数量化 II 類
11.数量化 III 類
12.数量化 IV 類
13.固有値・固有ベクトルを求める
14.逆行列を求める

---

堀　啓造(home page) e-mail hori@ec.kagawa-u.ac.jp