t分布の面積を計算する javascript(とガンマ関数)
カテゴリ:確率統計
function gamma(z, aqr){ var delta = 1.0/aqr; var area = 0; var oldArea = -1; var t = 0; while(t < 3 || Math.abs(area-oldArea) > Number.EPSILON){ oldArea = area; area += Math.pow(t, z-1)*Math.exp(-t)*delta; t += delta; } return area; } function tDist(df, z, aqr){ var halfN = (df+1)*0.5; var numerv = gamma(halfN, 1000); var denomv = gamma((df*0.5), 1000); var coeff = numerv/(Math.sqrt(Math.PI*df)*denomv); var invDF = 1/df; var delta = 1.0/aqr; var t = -z; var area = 0; var oldArea = -1; while(Math.abs(area-oldArea) > Number.EPSILON){ oldArea = area; area += coeff*Math.pow(1+t*t*invDF, -halfN)*delta; t += delta; } return area; }
引数の aqr は精度をきめる.ガンマ関数の aqr は 1000 ぐらいで十分だ. t分布の aqr は片側 0.0005 まで必要ならば,100,000 ぐらい必要になる. df は自由度で, tDist の z は t値だ.
しかしこの実装では実数の精度点で問題がある.
自由度1の面積が0になる.
精度 aqr が大きすぎると pow で精度が不足する. 自由度1の時だけ特別に計算式を別にする必要がある.
自由度が 300 以上の時,面積が無限大になる
ガンマ関数の pow で精度が不足するため t分布の係数を計算できない. そのため手動で係数を設定する必要がある.
上記の問題を対策したコード
function gamma(z, aqr){ var delta = 1.0/aqr; var area = 0; var oldArea = -1; var t = 0; if( z < 100){ while(t < 3 || Math.abs(area-oldArea) > Number.EPSILON){ oldArea = area; area += Math.pow(t, z-1)*Math.exp(-t)*delta; t += delta; } } else{ while(t < 3 || Math.abs(area-oldArea) > Number.EPSILON){ oldArea = area; var pPow = Math.pow(t, (z-1)*0.1); area += Math.exp(-t)*delta*pPow*pPow*pPow*pPow*pPow*pPow*pPow*pPow*pPow*pPow; t += delta; } } return area; } function tDist(df, z, aqr){ var halfN = (df+1)*0.5; var coeff = 0; if( df === 1 ){ coeff = 0.31831; } else if( 1 < df && df <= 300 ){ var numerv = gamma(halfN, 1000); var denomv = gamma((df*0.5), 1000); coeff = numerv/(Math.sqrt(Math.PI*df)*denomv); } else if( 300 < df && df <= 400 ){ coeff = 0.39863 + 0.00005*(df-300)/100; } else if( 400 < df && df <= 500 ){ coeff = 0.39868 + 0.00005*(df-400)/100; } else if( 500 < df && df < 1000000 ){ coeff = 0.39873 + 0.000212*(df-500)/1000000; } else{ coeff = 0.398942; } var invDF = 1/df; var delta = 1.0/aqr; var t = -z; var area = 0; var oldArea = -1; if( df === 1 ){ delta = 0.01; while(Math.abs(area-oldArea) > Number.EPSILON){ oldArea = area; area += coeff*delta/(1+t*t); t += delta; } } else{ while(Math.abs(area-oldArea) > Number.EPSILON){ oldArea = area; area += coeff*Math.pow(1+t*t*invDF, -halfN)*delta; t += delta; } } return area; }