use Math::Cephes qw(gamma);
	
# 非心カイ二乗分布 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] = 0.5 - $CumulativeDistributionFunction;
	# 上側累積確率 Upper Cumulative Distribution
	$NoncentralChiSquareDistribution[2] = 0.5 + $CumulativeDistributionFunction;
	
	return @NoncentralChiSquareDistribution;
}
	
# 確率密度関数 Frequency Function
# 引数 変数 自由度 非心度 ($X, $DegreesOfFreedom, $Noncentral)
# 戻り値 確率密度関数 ($FrequencyFunction)
sub FREQUENCYFUNCTION{
	my ($X, $DegreesOfFreedom, $Noncentral) = @_;
	my $FrequencyFunction = 0;
	my $Temp1 = (1 / 2) * exp(-(($X + $Noncentral) / 2));
	my $Temp2 = ($Noncentral != 0 ? ($X / $Noncentral) ** (($DegreesOfFreedom / 4) - (1 / 2)) : 0);
	my $BesselFunction = &BESSELFUNCTION((($DegreesOfFreedom / 2) - 1), sqrt($Noncentral * $X));
	
	# 非心カイ二乗分布 Frequency Function
	$FrequencyFunction = $Temp1 * $Temp2 * $BesselFunction;
	
	return $FrequencyFunction;
}
	
# 第１種変形ベッセル関数 Bessel Function
# 引数 次数 変数 ($N, $X)
# 戻り値 第１種変形ベッセル関数 ($BesselFunction)
sub BESSELFUNCTION{
	my ($N, $X) = @_;
	my $BesselFunction = 0;
	my $PrevBesselFunction = 0;
	my $Num1 = 0;
	my $Den1 = 0;
	my $Den2 = 0;
	my $Limit = 100;
	my $Epsilon = 1.0e-10;
	
	for(my $i = 0; $i < $Limit; $i++){
		# 分子
		$Num1 = ($X / 2) ** ($N + (2 * $i));
		# 分母
		$Den1 = ($i == 0 ? 1 : ($Den1 * $i));
		$Den2 = &gamma($N + $i + 1);
	
		# 一つ前
		$PrevBesselFunction = $BesselFunction;
		# 第１種変形ベッセル関数 Bessel Function
		$BesselFunction += ($Num1 / ($Den1 * $Den2));
	
		# 収束判定
		last if(abs($BesselFunction - $PrevBesselFunction) < $Epsilon);
	}
	
	return $BesselFunction;
}
	
# 近似 累積分布関数 Cumulative Distribution Function
# 引数 変数 自由度 非心度 ($X, $DegreesOfFreedom, $Noncentral)
# 戻り値 近似 累積分布関数 ($CumulativeDistributionFunction)
sub CUMULATIVEDISTRIBUTIONFUNCTION{
	my ($X, $DegreesOfFreedom, $Noncentral) = @_;
	my $CumulativeDistributionFunction = 0;
	my $D = $DegreesOfFreedom;
	my $N = $Noncentral;
	my $h = 1 - ( (2 / 3) * (( ($D + $N) * ($D + (3 * $N)) ) / ( ($D + (2 * $N)) * ($D + (2 * $N)) )));
	my $p = ( ($D + (2 * $N)) / (($D + $N) * ($D + $N)) );
	my $m = ($h - 1) * (1 - (3 * $h));
	my $Num1 = 1;
	my $Num2 = ($h * $p) * (1 - $h + (0.5 * (2 - $h) * ($m * $p)));
	my $Num3 = ($X / ($D + $N)) ** ($h);
	my $Den1 = $h * sqrt(2 * $p * (1 + ($m * $p)));
	
	# 近似 累積分布関数 Cumulative Distribution Function
	$CumulativeDistributionFunction = &SIMPSONSRULE(($Num1 - $Num2 - $Num3) / $Den1);
	
	return $CumulativeDistributionFunction;
}
	
# 確率密度関数 Frequency Function
# 引数 変数 ($X)
# 戻り値 確率密度関数 ($StandardNormalFrequencyFunction)
sub STANDARDNORMALFREQUENCYFUNCTION{
	my ($X) = @_;
	my $StandardNormalFrequencyFunction = 0;
	my $Pi = atan2(1, 1) * 4;
	
	# 標準正規分布 Frequency Function
	$StandardNormalFrequencyFunction = (1 / sqrt(2 * $Pi)) * exp(-(($X * $X) / 2));
	
	return $StandardNormalFrequencyFunction;
}
	
# シンプソンの公式 Simpson's Rule
# 引数 変数 ($X)
# 戻り値 シンプソンの公式 ($SimpsonsRule)
sub SIMPSONSRULE{
	my ($X) = @_;
	my $SimpsonsRule = 0;
	my $a = 0;
	my $b = $X;
	my $N = 100;
	my $h = ($b - $a) / $N;
	my $Sum = 0;
	my $SumA = 0;
	my $SumB = 0;
	
	# 計算
	$SumA = &STANDARDNORMALFREQUENCYFUNCTION($a);
	$SumB = &STANDARDNORMALFREQUENCYFUNCTION($b);
	
	# 分割数は適当
	for(my $i = 1; $i < $N; $i++){
		my $t = $a + ($i * $h);
		my $num = ($i % 2 == 0 ? 2 : 4);
		my $tmp = &STANDARDNORMALFREQUENCYFUNCTION($t);
	
		$Sum += $num * $tmp;
	}
	
	# シンプソンの公式 Simpson's Rule
	$SimpsonsRule = ($SumA + $Sum + $SumB) * ($h / 3);
	
	return $SimpsonsRule;
}