↑
広告
広告

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値だ．

しかしこの実装では実数の精度点で問題がある．

自由度１の面積が０になる．

精度 aqr が大きすぎると pow で精度が不足する． 自由度１の時だけ特別に計算式を別にする必要がある．

自由度が 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;
}

広告