ガンマ関数の実装とt分布

 「t分布」は統計で使用するものですが、統計学とはいかなるものか。

1.様々なデータを分析・比較し、傾向を見出すなどする学問。様々な学問や研究を下支えしている。
2.得体の知れない占いやインチキ科学などにおいて、それが科学的なものであるという印象を予備知識のない人々にアピールするためのキーワード。

 「砂山のパラドックス」という言葉があります。「砂山から1粒の砂を取り除いても、それは砂山だ。さらにもう1粒取り除いても同様だ。その調子で次々と砂山から砂を取り除き、砂山に1粒の砂しか残っていない状態になったら、それは砂山と呼べるのか」という問題です。
 この問題に対する答えは、おそらく多種多様なものになるはずです。「1粒でも砂がある以上は砂山と定義できる」と考える人がいれば、「n粒以上の砂があるのが砂山だ」と条件を定める人、「見た目が砂山であるものが砂山だ」と考える人もいるでしょう。それでは結局、何が砂山で何がそうではないのでしょうか。
 さて、世の中は嫌煙ムードで、タバコが体に毒であるのはもはや常識となっています。ところが、未だにそれを認めようとしない人もごくごく少数ながら存在するようです。そのような人が自説を声高に主張するせいで、非喫煙者を怒らせて余計に喫煙者の肩身を狭くしている部分があり、他人の迷惑を慮れるまっとうな喫煙者の方はさぞ迷惑をしておられるのではないかと考えたりします。
 仮にここに「喫煙者は肺ガンが2倍」という統計データがあるとしましょう。ところが、「嫌煙ファシズム」などと意味の分からないことを称している人と討論しようものなら、以下のような意味不明な経緯をたどりかねません。

「この通り、喫煙者は肺ガンが倍増します。これでも健康に害はないとおっしゃいますか?」
「サンプル数はいくつだ。少ないのではないか」
「サンプルは100万ですから、決して少ない数ではないでしょう」
「たまたま喫煙者のサンプルに肺ガンの人が多かったのだろう」
「100万も調査して、肺ガンが2倍であっても、たまたまなのですか?」
「当たり前だ。どうして偶然2倍にならないと言い切れる?たまたまではない証拠があるか?」
「それでは、日本人全員を調査しろとおっしゃるのですか?」
「日本人全員?バカ言っちゃいけない。今の日本人全員だって、たまたま喫煙者に肺ガンが多いかもしれないじゃないか」
「ではどうすればよろしいのですか。逆に、喫煙が肺ガンを生じさせない証拠でもあるのですか?」
「悪魔の証明だ!証拠がないなら、害はないとするのが当たり前だろう!」

 この種の自分勝手な人は、たとえ病気の発生率が500倍であったとしても、決して害を認めようとはしないでしょう。「500倍も偶然のうち」なのですから。
 それでは、実際問題として一体何倍までであれば偶然で、何倍以上なら偶然ではないと定義できるのでしょうか。これを知るための方法として、統計学には「信頼区間」というものがあり、限られたサンプルのデータから全体の平均値を予想することができます。例えば「喫煙者全員の発ガン率の平均は、非喫煙者の1.95〜2.05倍の範囲内にある確率が95%である」「1.85〜2.15倍の範囲にある確率が99%である」といった具合です。もし「それでは、その区間に含まれない可能性もあるわけだ」と反論を受けたとしても、外れる可能性はせいぜい1%や5%ですから、「発ガン性がある可能性が極めて高いという結果が出た以上、これが誤りと言うなら発ガン性がない証明をしてください」と返せるわけです。なお、ここでは喫煙を例に出しましたが、あくまで身近な例として用いただけであって、別にタバコの害について論じる気はないこと、値は全く根拠のない架空のものであることをお断りしておきます。
 統計学は科学その他の証明の他、疾患に対して新治療が有意に効果をもたらすかを調べたり、安全性が要求される機械の耐久性を調べたり、社会問題を統計によって明らかにしたりと、実に様々な分野で活用されています。欠かせない学問といえるでしょう。
 ところが、これだけ価値ある学問とあれば、ネームバリューにただ乗りしようとする人間は出てくるものです。占いや迷信、疑似科学の紹介に踊る「統計学」の文字の何とむなしいことか。統計学を名乗るからには、

・使用したデータの概要。サンプルは十分か、回答者が偏る要素はないか、設問は公平か。「占いの本を買った人が返送したアンケートの集計」を根拠に血液型占いの正当性を訴える人もいるようですが、その集計には何の価値もありません。「おもいッきり」や「あるある」などに至っては、サンプルが2つや4つでは全く話になりません。
・分析方法の開示。同じ条件での追試で同等の結果が出ることが重要です。分析方法を開示せず、単に「統計学」を標榜してもっともらしいことを言う主張は、統計学の名をかたることで信憑性を上げようとした疑似科学であって、統計学の冒涜に他なりません。

 この程度は必要でしょう。逆に、これらが欠如しているばかりか、単に「統計学」としか書いていない場合は、その迷信は統計学の名を不当に利用していると考えてほぼ間違いありません。学問たるものを何だと考えているのか、実に許しがたいことです。阪神優勝を毎年予想していたらそのうち見事的中した、とある日本一有名なペテン師の占いもこれに含まれます。

 統計学で有名なものに、スチューデントのt分布やt検定というものがあります。tならティーチャーのt検定でもよさそうなものですが、発表者のペンネームから来ているようです。何やら難解そうですが、使うだけなら誰にでもできます。表計算ソフトにはTDISTやTTEST、TINVといった関数が用意されており、統計を簡単に試してみるにはこの上ないツールです。
 以下、簡単なサンプルを。

問題
 とある機械装置について、10個を無作為に抽出して耐久性の調査を行いました。その結果は以下の通りです。

1015 , 988 , 992 , 975 , 1022 , 1001 , 998 , 1008 , 1011 , 994

  1. サンプル10個の耐久性平均値は1000.4ですが、機械装置全体の耐久性の平均値はサンプルの平均値からどの程度の範囲内に位置すると推測されますか。t分布を用い、信頼区間95%及び99%の場合の結果を示してください。
  2. もし機械装置全体の耐久性の平均が975以下になってしまった場合、この装置は赤字製品となることが分かっています。高い確率で975以下となるようなら、発売を見送らなくてはなりません。平均値が975以下となる確率はどの程度ですか。t分布を用いて調べてください。
答え
 おそらく次のようになります。

- A B
1Sample 11015
2Sample 2988
3Sample 3992
4Sample 4975
5Sample 51022
6Sample 61001
7Sample 7998
8Sample 81008
9Sample 91011
10Sample 10994
11 平均 =AVERAGE(B1:B10)
12 95%範囲 =TINV(0.05,COUNT(B$1:B$10)-1)*SQRT(VAR(B$1:B$10)/COUNT(B$1:B$10))
13 99%範囲 =TINV(0.01,COUNT(B$1:B$10)-1)*SQRT(VAR(B$1:B$10)/COUNT(B$1:B$10))
14 975以下の確率 =TDIST(ABS(B$11-975)/SQRT(VAR(B$1:B$10)/COUNT(B$1:B$10)),COUNT(B$1:B$10)-1,1)

 TDISTがt分布で、TINVが逆関数です。SQRT(VAR(...)/COUNT(...))は二乗平均平方根で、解説サイトによっては「標準偏差」と書いてあったりしますが、ExcelのSTDEVやSQLのSTDDEVとは別物です。TDISTの3つ目の引数は1か2しか指定できず、今回の状況で2にすると「全体の平均値がABS(1000.4-975)、つまり+-25.4以上離れる確率」となってしまいます。今回は「1000.4 - 25.4」に限定したいので、1を使っています。
 私の環境では、結果は次のようになりました。

- A B
11 平均 1000.4
12 95%範囲 10.02181
13 99%範囲 14.39745
14 975以下の確率 0.000141

 つまり、機械装置の耐久性の平均は95%の確率で1000.4 +-10.02181の範囲内となり、99%の確率で+-14.39745の範囲内となります。また、耐久性の平均が975以下となってしまう確率は0.014%程度ですから、発売してもほぼ問題ないといえます(現実問題として、サンプルが10個では「あるある」並みに信用できないのですが、ここでは問いません)。

 Excelでt分布の関数があるのなら、他の言語でも作れるはずです。Wikipediaで算出方法を調べてみると、t分布の項にはガンマ関数が必要であるとあります。したがって、t分布の実装を作るには、ガンマ関数を実装しなくてはならないようです。
 それではガンマ関数とは何か。「5!」や「10!」は「階乗」と呼ばれますが、ガンマ関数は「3.5」や「0.1」などにも階乗を適用できるようにするものです。ただし階乗とは数字が1つずれて、「gamma(n) = (n - 1)!」になるとのこと。定義は上記Wikipediaの項目にありますが、関数言語風に言えば以下のようになります。
let gamma z = integral 0 infinity (fun t -> (t ** (z - 1)) * (exp(1) ** -t))
(* integral は
(start : float) (end : float) (d : float -> float)
のような架空の関数とする *)
 式の見た目は単純ですが、不可解な曲線を描く積分というものはとにかく厄介で、しかも無限大まで続くのですからたまりません。無限の細かさを持つものを、さらに無限大に至るまで計算しているわけにはいきませんし、かといって制度は保たなくてはなりません。
 日本語Wikipediaにはスターリングの近似なる項目があり、ここにある式は比較的簡単ながらそこそこの近似値が得られるようになっています。ただ、近似値の精度は(方法にもよるとはいえ)Excelのgammalnに劣り、納得できる水準ではありません。
 そこで使うのがLanczos approximation(ランチョス近似)なのですが、何と日本語Wikipediaには今のところ記事がありません。英語版Wikipediaの項目にしても、実装するには何をどうしたものかさっぱり分かりません。
 しかしながら、その外部リンクにはThe Lanczos Approximationという非常に優良な資料があります。なぜかマトリックスが出てきており、これがガンマ関数と何の関係があるのか私にはさっぱり理解できませんが、とりあえず書いてあることをすべてPHPでそのまま実装してみたのが以下のプログラムです。
<?php
function fact($n){
	$result = 1;
	for($i = 2; $i <= $n; $i++)
		$result *= $i;

	return $result;
}

function combination($n , $k){
	if($k < 0)
		return 0;
	return fact($n) / (fact($k) * fact($n - $k));
}

class Matrix{
	private $matrix;
	private $height;
	private $width;

	function __construct($y , $x){
		$this->height = $y;
		$this->width = $x;
		$this->matrix = array_fill(0 , $y , array_fill(0 , $x , 0));
	}

	static function createNeutral($length){
		$matrix = new Matrix($length , $length);
		for($y = 0; $y < $length; $y++){
			for($x = 0; $x < $length; $x++)
				$matrix->set($y , $x , $y == $x ? 1 : 0);
		}

		return $matrix;
	}
	static function createVector($source){
		$length = count($source);
		$matrix = new Matrix(1 , $length);
		for($x = 0; $x < $length; $x++)
			$matrix->set(0 , $x , $source[$x]);

		return $matrix;
	}
	function rotate(){
		$m = new Matrix($this->getWidth() , $this->getHeight());
		for($y = 0; $y < $this->getHeight(); $y++){
			for($x = 0; $x < $this->getWidth(); $x++)
				$m->set($x , $y , $this->get($y , $x));
		}

		return $m;
	}

	function get($y , $x){
		if(($this->matrix[$y][$x]) === null)
			print "NULL!";
		return $this->matrix[$y][$x];
	}
	function set($y , $x , $v){
		$this->matrix[$y][$x] = $v;
	}
	function getWidth(){
		return $this->width;
	}
	function getHeight(){
		return $this->height;
	}

	function multiply($n){
		$m = $this;

		if($m->getWidth() != $n->getHeight())
			return null;

		$width = $n->getWidth();
		$height = $m->getHeight();

		$matrix = new Matrix($height , $width);

		for($ay = 0; $ay < $height; $ay++){
			for($ax = 0; $ax < $width; $ax++){
				// 演算
				$value = 0;
				for($i = 0; $i < $m->getWidth(); $i++)
					$value += $m->get($ay , $i) *
						$n->get($i , $ax);
				$matrix->set($ay , $ax , $value);
			}
		}

		return $matrix;
	}
}

function show_matrix(Matrix $matrix){
	print "<table border="1">n";
	for($y = 0; $y < $matrix->getHeight(); $y++){
		print "<tr>";
		for($x = 0; $x < $matrix->getWidth(); $x++){
			print "<td>";
			print $matrix->get($y , $x);
			print "</td>";
		}
		print "</tr>n";
	}
	print "</table>";
}

function gammaln($z , $n = 8 , $g = null){
	if($g == null)
		$g = $n - 0.5;

	$z --;

	// B 作成
	$B = new Matrix($n , $n);
	for($i = 0; $i < $n; $i++){
		for($j = 0; $j < $n; $j++){
			$value = 0;
			if($i == 0)
				$value = 1;
			else if($i > 0 && j >= i)
				$value = pow(-1 , $j - $i) *
					combination($i + $j - 1 , $j - $i);
			else
				$value = 0;
			$B->set($i , $j , $value);
		}
	}

	// C 作成
	$C = new Matrix($n , $n);
	for($i = 0; $i < $n; $i++){
		for($j = 0; $j < $n; $j++){
			$value = 0;
			if($i == 0 && $j == 0){
				$value = 0.5;
			}else if($j > $i){
				$value = 0;
			}else{
				$sigma = 0;
				for($k = 0; $k <= $i; $k++){
					$sigma += combination(2 * $i , 2 * $k) *
						combination($k , $k + $j - $i);
				}
				$value = pow(-1 , $i - $j) * $sigma;
			}
			$C->set($i , $j , $value);
		}
	}

	// D 作成
	$D = new Matrix($n , $n);
	for($i = 0; $i < $n; $i++){
		for($j = 0; $j < $n; $j++){
			$value = 0;
			if($i != $j)
				$value = 0;
			else if($i == 0 && $j == 0)
				$value = 1;
			else if($i == 1 && $j == 1)
				$value = -1;
			else
				$value = ($D->get($i - 1 , $j - 1) *
					(2 * (2 * $i - 1))) / ($i - 1);
			$D->set($i , $j , $value);
		}
	}

	// F 作成
	$F = new Matrix($n , 1);
	for($i = 0; $i < $n; $i++){
		$value = (fact(2 * $i) * exp($i + $g + 0.5)) /
			(fact($i) * pow(2 , 2 * $i - 1) *
			pow($i + $g + 0.5 , $i + 0.5));
		$F->set($i , 0 , $value);
	}

	// P 作成
	$P = $D->multiply($B)->multiply($C)->multiply($F);

	// Z 作成
	$Z = new Matrix(1 , $n);
	$Z->set(0 , 0 , 1);
	for($i = 1; $i < $n; $i++)
		$Z->set(0 , $i , 1 / ($z + $i));

	print "D = ";
	show_matrix($D);

	print "B = ";
	show_matrix($B);

	print "C = ";
	show_matrix($C);

	print "F = ";
	show_matrix($F);

	print "Z = ";
	show_matrix($Z);

	print "P = ";
	show_matrix($P);

	$r = log($Z->multiply($P)->get(0 , 0)) + ($z + 0.5) *
		log($z + $g + 0.5) - ($z + $g + 0.5);
	print "gammaln = " . $r;
	return $r;
}

// 使用例
print " , gamma = " . exp(gammaln(-3/2));
?>
 上記gammaln関数が取る$nは「精度」のようなもので、何らかの整数を指定します。$gは「任意の値」とされているようで、好きな実数を与えて構わないようですが、どこかの資料で「gはn - 0.5程度が無難」との記述を読んだように記憶しているため、実装では指定のない限りそうしています。実際にはnの値によって最適なgの値があるようです。
 これを実行してみると、以下のような結果が得られます。n = 8 , g = 7.5です。
D =
10000000
0-1000000
00-600000
000-300000
0000-140000
00000-63000
000000-27720
0000000-12012
B =
11111111
-01-23-45-67
0-01-410-2035-56
-00-01-621-56126
0-00-01-836-120
-00-00-01-1055
0-00-00-01-12
-00-00-00-01
C =
0.50000000
-12000000
1-8800000
-118-48320000
1-32160-256128000
-150-4001120-128051200
1-72840-35846912-614420480
-198-15689408-2688039424-286728192
F =
2107.85560707
300.11421954
104.480701073
50.8623205503
29.7383924002
19.5187343297
13.8664417012
10.4308729974
Z =
1-0.666666666667-220.6666666666670.40.2857142857140.222222222222
P =
2.50662827329
2901.41861668
-5929.16815028
4148.27391368
-1164.7605648
117.413508415
-2.78664112091
0.00379359722137
gammaln = 0.860047015803 , gamma = 2.36327180221
 Wikipediaを見てみても、ガンマ「-2/3」の近似値は2.363とありますので、正しく動作しているようです。また、整数を与えてみても階乗と同じ結果(引数の値は1ずつずれるので注意)が得られます。
 ここで面白いのは、gammaln内でいくつものマトリックスを作成していますが、$Z以外のマトリックスにはgammalnの引数zが必要ない点です。gammalnの引数n及びgは精度及び任意の定数ですので、引数zが1でも10でも-0.5でも関係ありません。したがって、精度を可変にする必要がなければ、$Pは定数にできます。
 $Pも$Zもマトリックスの扱いではありますが、実際にはどちらもベクトルとなっていますので、掛け算は単純にそれぞれの各数値を乗算するだけで構いません。つまり、引数から$Zを求めて$Pの定数と掛け合わせてしまえば、それでgammaln関数は実装できてしまいます。
 以下、gammaln(z , 8 , 7.5)を用いた場合の$Pを定数として使ったgammaln関数の実装です。
<?php
header("Content-type: text/plain");

function gammaln($z){
	$z --;

	$p = array(2.50662827329 , 2901.41861668 , -5929.16815028 , 4148.27391368 ,
		-1164.7605648 , 117.413508415 , -2.78664112091 , 0.00379359722137);
	$g = 7.5;

	$result = $p[0];
	for($i = 1; $i < count($p); $i++)
		$result += $p[$i] * (1 / ($z + $i));

	return log($result) + ($z + 0.5) * log($z + $g + 0.5) - ($z + $g + 0.5);
}

print exp(gammaln(-3/2));
?>
 何とも短いコードとなってしまいました。
 ガンマ関数さえ実装できれば、t分布を求める関数は簡単に作成できます。
<?php
function gammaln($z){
	$z --;

	$p = array(2.50662827329 , 2901.41861668 , -5929.16815028 , 4148.27391368 ,
		-1164.7605648 , 117.413508415 , -2.78664112091 , 0.00379359722137);
	$g = 7.5;

	$result = $p[0];
	for($i = 1; $i < count($p); $i++)
		$result += $p[$i] * (1 / ($z + $i));

	return log($result) + ($z + 0.5) * log($z + $g + 0.5) - ($z + $g + 0.5);
}

function gamma($z){
	return exp(gammaln($z));
}

// 補足:ベータ関数の実装
function beta($x , $y){
	return gamma($x) * gamma($y) / gamma($x + $y);
}

// t分布のPDF
function tpdf($n , $t){
	$f = $n - 1;
	$a = gamma(($f + 1) / 2);
	$b = sqrt($f * pi()) * gamma($f / 2);
	$c = pow(1 + pow($t , 2) / $f , -(($f + 1) / 2));
	return $a / $b * $c;

	// ベータ関数を用い、以下のようにしても同じ結果が得られる
	//return (1 / (sqrt($v) * beta(0.5 , $v / 2))) * pow(1 + $t * $t / $v , -(($v + 1) / 2));
	// あるいは
	//return pow($v / ($v + $t * $t) , (1 + $v) / 2) / (sqrt($v) * beta($v / 2 , 0.5));
}

// gd2を使用してグラフを作成してみる
$width = 400;
$height = 400;

$img = imagecreatetruecolor($width , $height);
imagefilledrectangle($img , 0 , 0 , $width - 1 , $height - 1 , 255 << 16 | 255 << 8 | 255);

$area = 4;

$frees = array(2 , 5 , 10 , 30 , 100 , 300);
$colors = array(255 << 16 , (255 << 16) | 255 , (255 << 16) | (255 << 8) , 255 << 8 ,
	(255 << 8) | 255 , 255);

for($a = 0; $a < count($frees); $a++){
	$f = $frees[$a];

	$before = tpdf($f , -$area) * 1000;
	for($i = 1; $i < $width; $i++){
		$t = $area * 2 * $i / $width - $area;
		$result = tpdf($f , $t) * 1000;

		imageline($img , $i - 1 , $height - $before , $i , $height - $result , $colors[$a]);

		$before = $result;
	}
}

header("Content-type: image/gif");
imagegif($img);
imagedestroy($img);
?>
 目盛りなどがないため少々見づらいですが、これでt分布のグラフが出力されます。赤が自由度2、青が自由度300の場合のt分布です。



 自由度が高いほど中央が大きく、そうでないほど裾野が広くなります。