# スネデカーのF分布 Snedecor's F-Distribution
# 引数 変数 分子の自由度 分母の自由度 ($X, $NumDegreesOfFreedom, $DenDegreesOfFreedom)
# 戻り値 スチューデントのt分布 (@SnedecorsFDistribution)
sub SNEDECORSFDISTRIBUTION{
	my ($X, $NumDegreesOfFreedom, $DenDegreesOfFreedom) = @_;
	my @SnedecorsFDistribution = ();
	my $IncompleteBetaFunction = 0;
	my $t = 0;
	
	# 変数,自由度の確認
	if(($X < 0) || ($NumDegreesOfFreedom <= 0) || ($DenDegreesOfFreedom <= 0)){
		return 0;
	}
	
	# 確率密度関数 Frequency Function
	$SnedecorsFDistribution[0] = &FREQUENCYFUNCTION($X, $NumDegreesOfFreedom, $DenDegreesOfFreedom);
	
	$t = (($NumDegreesOfFreedom * $X) / (($NumDegreesOfFreedom * $X) + $DenDegreesOfFreedom));
	# 累積分布関数 正規化不完全ベータ関数 Incomplete Beta Function
	$IncompleteBetaFunction = &REGULARIZEDINCOMPLETEBETAFUNCTION(($NumDegreesOfFreedom / 2), ($DenDegreesOfFreedom / 2), $t);
	
	# 下側累積確率 Lower Cumulative Distribution
	$SnedecorsFDistribution[1] = $IncompleteBetaFunction;
	# 上側累積確率 Upper Cumulative Distribution
	$SnedecorsFDistribution[2] = 1 - $IncompleteBetaFunction;
	
	return @SnedecorsFDistribution;
}
	
# 確率密度関数 Frequency Function
# 引数 変数 分子の自由度 分母の自由度 ($X, $NumDegreesOfFreedom, $DenDegreesOfFreedom)
# 戻り値 確率密度関数 ($FrequencyFunction)
sub FREQUENCYFUNCTION{
	my ($X, $NumDegreesOfFreedom, $DenDegreesOfFreedom) = @_;
	my $FrequencyFunction = 0;
	my $d1 = $NumDegreesOfFreedom;
	my $d2 = $DenDegreesOfFreedom;
	my $Num1 = 0;
	my $Num2 = 0;
	my $Num3 = 0;
	my $Den = 0;
	
	# 分子
	$Num1 = (($d1 * $X) / (($d1 * $X) + $d2)) ** ($d1 / 2);
	$Num2 = (1 - (($d1 * $X) / (($d1 * $X) + $d2))) ** ($d2 / 2);
	$Num3 = $X ** (-1);
	# 分母
	$Den = &BETAFUNCTION(($d1 / 2), ($d2 / 2), 1);
	# スネデカーのF分布 Snedecor's F-Distribution Frequency Function
	$FrequencyFunction = ($Num1 * $Num2 * $Num3) / $Den;
	
	return $FrequencyFunction;
}

# 正則化された不完全ベータ関数 Regularized Incomplete Beta Function
# 引数 変数 変数 変数 ($AA, $BB, $X)
# 戻り値 正則化された不完全ベータ関数 ($RegularizedRegularizedIncompleteBetaFunction)
sub REGULARIZEDINCOMPLETEBETAFUNCTION{
	my ($AA, $BB, $X) = @_;
	my $BetaFunction = 0;
	my $IncompleteBetaFunction = 0;
	my $RegularizedIncompleteBetaFunction = 0;

	# ベータ関数 Beta Function
	$BetaFunction = &BETAFUNCTION($AA, $BB, 1);
	# 不完全ベータ関数 Incomplete Beta Function
	$IncompleteBetaFunction = &BETAFUNCTION($AA, $BB, $X);

	# 正則化された不完全ベータ関数 Regularized Incomplete Beta Function
	$RegularizedIncompleteBetaFunction = $IncompleteBetaFunction / $BetaFunction;

	return $RegularizedIncompleteBetaFunction;
}

# ベータ関数 Beta Function
# 引数 変数 変数 変数 ($AA, $BB, $X)
# 戻り値 ベータ関数 (@BetaFunction)
sub BETAFUNCTION{
	my ($AA, $BB, $X) = @_;
	my $BetaFunction = 0;
	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;
	
	# シンプソンの公式 Simpson's Rule
	$SumA = ($a ** ($AA - 1)) * ((1 - $a) ** ($BB - 1));
	$SumB = ($b ** ($AA - 1)) * ((1 - $b) ** ($BB - 1));
	
	# 分割数は適当
	for(my $i = 1; $i < $N; $i++){
		my $t = $a + ($i * $h);
		my $num = ($i % 2 == 0 ? 2 : 4);
		my $tmp = ($t ** ($AA - 1)) * ((1 - $t) ** ($BB - 1));
	
		$Sum += $num * $tmp;
	}
	
	# ベータ関数 Beta Function
	$BetaFunction = ($SumA + $Sum + $SumB) * ($h / 3);
	
	return $BetaFunction;
}