FC2ブログ
最初
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/NoncentralChiSquareDistribution1.txt

修正1 8/20

# 非心カイ二乗分布 Noncentral Chi-Square Distribution
# 引数 変数 自由度 非心度 ($X, $DegreesOfFreedom, $Noncentral)
# 戻り値 非心カイ二乗分布 (@NoncentralChiSquareDistribution)
sub NONCENTRALCHISQUAREDISTRIBUTION{
my ($X, $DegreesOfFreedom, $Noncentral) = @_;
my @NoncentralChiSquareDistribution = ();
my $CumulativeDistributionFunction = 0;

# 変数,自由度,非心度の確認
if(($X < 0) || ($DegreesOfFreedom <= 0) || ($Noncentral < 0)){
return 0;
}

# 確率密度関数 Frequency Function
$NoncentralChiSquareDistribution[0] = &FREQUENCYFUNCTION($X, $DegreesOfFreedom, $Noncentral);

# 近似 累積分布関数 Cumulative Distribution Function
$CumulativeDistributionFunction = &CUMULATIVEDISTRIBUTIONFUNCTION($X, $DegreesOfFreedom, $Noncentral);

# 下側累積確率 Lower Cumulative Distribution
$NoncentralChiSquareDistribution[1] = $CumulativeDistributionFunction;
# 上側累積確率 Upper Cumulative Distribution
$NoncentralChiSquareDistribution[2] = 1 - $CumulativeDistributionFunction;

return @NoncentralChiSquareDistribution;
}

# 確率密度関数 Frequency Function
# 引数 変数 自由度 非心度 ($X, $DegreesOfFreedom, $Noncentral)
# 戻り値 確率密度関数 ($FrequencyFunction)
sub FREQUENCYFUNCTION{
my ($X, $DegreesOfFreedom, $Noncentral) = @_;
my $FrequencyFunction = 0;
my $Sum = 0;
my $PrevSum = 0;
my $Num1 = 0;
my $Num2 = 0;
my $Num3 = 0;
my $Den1 = 1;
my $Den2 = 0;
my $Den3 = 0;
my $Temp = 0;
my $Limit = 100;
my $Epsilon = 1.0e-10;

for(my $i = 0; $i < $Limit; $i++){
$Temp = ($DegreesOfFreedom + (2 * $i)) / 2;

# 分子
$Num1 = exp(-($Noncentral / 2)) * (($Noncentral / 2) ** $i);
$Num2 = $X ** ($Temp - 1);
$Num3 = exp(-($X / 2));
# 分母
$Den1 *= $i if($i != 0);
$Den2 = 2 ** $Temp;
$Den3 = &GAMMAFUNCTION($Temp);

$PrevSum = $Sum;
$Sum += (($Num1 * $Num2 * $Num3) / ($Den1 * $Den2 * $Den3));

# 収束判定
last if(abs($Sum - $PrevSum) < $Epsilon);
}

# 非心カイ二乗分布 Frequency Function
$FrequencyFunction = $Sum;

return $FrequencyFunction;
}

# 累積分布関数 Cumulative Distribution Function
# 引数 変数 自由度 非心度 ($X, $DegreesOfFreedom, $Noncentral)
# 戻り値 累積分布関数 ($CumulativeDistributionFunction)
sub CUMULATIVEDISTRIBUTIONFUNCTION{
my ($X, $DegreesOfFreedom, $Noncentral) = @_;
my $CumulativeDistributionFunction = 0;
my $Sum = 0;
my $PrevSum = 0;
my $Temp1 = 0;
my $Temp2 = 0;
my $Temp3 = 0;
my $Factorial = 1;
my $Limit = 100;
my $Epsilon = 1.0e-10;

for(my $i = 0; $i < $Limit; $i++){
$Factorial *= $i if($i != 0);
$Temp1 = exp(-($Noncentral / 2));
$Temp2 = (($Noncentral / 2) ** $i) / $Factorial;
$Temp3 = &COMPLEMENTGAMMAFUNCTIONFIRST((($DegreesOfFreedom + (2 * $i)) / 2),($X / 2)) / &GAMMAFUNCTION(($DegreesOfFreedom + (2 * $i)) / 2);

$PrevSum = $Sum;
$Sum += ($Temp1 * $Temp2 * $Temp3);

# 収束判定
last if(abs($Sum - $PrevSum) < $Epsilon);
}

# 累積分布関数 Cumulative Distribution Function
$CumulativeDistributionFunction = $Sum;

return $CumulativeDistributionFunction;
}

# ガンマ関数(近似) Gamma Function Approximation
# 引数 値 ($X)
# 戻り値 ガンマ関数(近似) ($GammaFunction)
sub GAMMAFUNCTION{
my ($X) = @_;
my $GammaFunction = 0;
my $Factorial = 1;
my $Diff = abs($X - int($X));
my $Pi = atan2(1, 1) * 4;
my $Temp = 0;

if($Diff == 0){
for(my $i = $X - 1; $i >= 2; $i--){
$Factorial *= $i;
}

# X-1の階乗 ガンマ関数 Gamma Function
$GammaFunction = $Factorial;
}elsif($Diff == 0.5){
my $Start = ($X > 0 ? ((2 * int($X)) - 1) : ((2 * int(abs($X) + 1)) - 1));

for(my $i = $Start; $i > 1; $i -= 2){
$Factorial *= $i;
}

if($X > 0){
$Temp = ($Factorial / (2 ** int($X))) * sqrt($Pi);
}else {
$Temp = ((((-1) ** int(abs($X) + 1)) * (2 ** int(abs($X) + 1))) / $Factorial) * sqrt($Pi);
}

# 半整数 ガンマ関数 Gamma Function
$GammaFunction = $Temp
}else {
if($X > 0){
$Temp = exp(&LOGGAMMAFUNCTION($X));
}else {
$Temp = $Pi / (sin($Pi * $X) * exp(&LOGGAMMAFUNCTION(1 - $X)));
}

# ガンマ関数(近似) Gamma Function Approximation
$GammaFunction = $Temp
}

return $GammaFunction;
}

# ログガンマ関数(近似) Log Gamma Function Approximation
# 引数 値 ($x)
# 戻り値 ログガンマ関数(近似) ($LogGammaFunction)
sub LOGGAMMAFUNCTION{
my ($x) = @_;
my $LogGammaFunction = 0;
my @Bernoulli = ((1 / 12), (1 / 360), (1 / 1260), (1 / 1680), (1 / 1188), (691 / 360360), (1 / 156), (3617 / 122400), (43867 / 244188), (174611 / 125400), (77683 / 5796), (236364091 / 1506960));
my $Pi = atan2(1, 1) * 4;
my $X = $x;
my $Sum = 0;
my $Sign = 0;
my $Power1 = 1;
my $Power2 = 1;
my $Count = @Bernoulli - 1;

for($X = $x; $X <= $Count; $X++){
$Power1 *= $X;
}
$Power2 = 1 / ($X * $X);

for(my $i = $Count; $i >= 0; $i--){
# 符号
$Sign = (($i % 2) == 0 ? 1: -1);

$Sum = $Sum + ($Sign * $Bernoulli[$i]);
$Sum = $Sum * $Power2 if($i != 0);
}
$Sum = $Sum / $X;

# ログガンマ関数(近似) Log Gamma Function Approximation
$LogGammaFunction = ((1 / 2) * log(2 * $Pi)) - log($Power1) - $X + (($X - (1 / 2)) * log($X)) + $Sum;

return $LogGammaFunction;
}

# 第1種不完全ガンマ関数 Incomplete Gamma Function Of The First Kind
# 引数 値 値 ($A, $X)
# 戻り値 第1種不完全ガンマ関数 ($IncompleteGammaFunction)
sub COMPLEMENTGAMMAFUNCTIONFIRST{
my ($A, $X) = @_;
my $IncompleteGammaFunction = 0;
my $IncompleteGamma = 0;
my $PrevIncompleteGamma = 0;
my $Factorial = 1;
my $Num = 0;
my $Den = 1;
my $Sum = 0;
my $Sign = 0;
my $Limit = 100;
my $Epsilon = 1.0e-10;

# 値の確認
if($X < 0){
return 0;
}

if(($A > 0) && (($A - int($A)) == 0)){
# 正の整数
for(my $i = 0; $i < $A; $i++){
# 階乗
$Factorial = ($i == 0 ? 1 : $Factorial * $i);
# 分子
$Num = $X ** $i;
# 分母
$Den = $Factorial;

$Sum += ($Num / $Den);
}

# 第1種不完全ガンマ関数 Incomplete Gamma Function Of The First Kind
$IncompleteGamma = $Factorial * (1 - (exp(-$X) * $Sum));
}else {
for(my $i = 0; $i < $Limit; $i++){
# 符号
$Sign = (($i % 2) == 0 ? 1: -1);
# 分子
$Num = $X ** $i;
# 分母
$Den = $A + $i;
# 階乗
$Factorial = ($i == 0 ? 1 : $Factorial * $i);

# 一つ前
$PrevIncompleteGamma = $IncompleteGamma;
# 第1種不完全ガンマ関数 Incomplete Gamma Function Of The First Kind
$IncompleteGamma += ($Sign * ($Num / ($Den * $Factorial)));

# 収束判定
last if(abs($IncompleteGamma - $PrevIncompleteGamma) < $Epsilon);
}

# 第1種不完全ガンマ関数 Incomplete Gamma Function Of The First Kind
$IncompleteGamma = ($X ** $A) * $IncompleteGamma
}

# 第1種不完全ガンマ関数 Incomplete Gamma Function Of The First Kind
$IncompleteGammaFunction = $IncompleteGamma;

return $IncompleteGammaFunction;
}


第1種ベッセル関数
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/BesselFunction1.txt
第1種変形ベッセル関数
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/BesselFunction2.txt
第1種ベッセル関数(漸化式)
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/BesselFunction3.txt

参考URL
Noncentral chi-square distribution - Wikipedia, the free encyclopedia
非心カイ2乗分布 - 高精度計算サイト
第1種ベッセル関数(グラフ) - 高精度計算サイト
第1種変形ベッセル関数(グラフ) - 高精度計算サイト
Bessel Function

一言
wikipediaに載っている累積分布関数では上手く値が出せなかったので、近似を使用
非心度を0にしカイ二乗分布を出すと確率密度が0になる
第1種変形ベッセル関数の漸化式は分からん
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

Author:雨宮
Firefoxを使用しているので気づかなかったけど、IE6でソースコードを上手くコピーできない

5/3
携帯用ならIE6でもソースコードをコピーできる
携帯用

検索フォーム


あわせて読みたいブログパーツ
一寸先は闇 RSS