# スチューデントのt分布 Student's T-Distribution
# 引数 変数 自由度 ($X, $DegreesOfFreedom)
# 戻り値 スチューデントのt分布 (@StudentsTDistribution)
sub STUDENTSTDISTRIBUTION{
	my ($X, $DegreesOfFreedom) = @_;
	my @StudentsTDistribution = ();
	my $IncompleteBetaFunction = 0;
	my $t = 0;
	
	# 自由度の確認
	if($DegreesOfFreedom <= 0){
		return 0;
	}
	
	# 確率密度関数 Frequency Function
	$StudentsTDistribution[0] = &FREQUENCYFUNCTION($X, $DegreesOfFreedom);
	
	$t = ($X + sqrt(($X * $X) + $DegreesOfFreedom)) / (2 * sqrt(($X * $X) + $DegreesOfFreedom));
	# 累積分布関数 正則不完全ベータ関数 Regularized Incomplete Beta Function
	$IncompleteBetaFunction = &REGULARIZEDINCOMPLETEBETAFUNCTION(($DegreesOfFreedom / 2), ($DegreesOfFreedom / 2), $t);
	
	# 下側累積確率 Lower Cumulative Distribution
	$StudentsTDistribution[1] = $IncompleteBetaFunction;
	# 上側累積確率 Upper Cumulative Distribution
	$StudentsTDistribution[2] = 1 - $IncompleteBetaFunction;
	
	return @StudentsTDistribution;
}
	
# 確率密度関数 Frequency Function
# 引数 変数 自由度 ($X, $DegreesOfFreedom)
# 戻り値 確率密度関数 ($FrequencyFunction)
sub FREQUENCYFUNCTION{
	my ($X, $DegreesOfFreedom) = @_;
	my $FrequencyFunction = 0;
	my $Pi = atan2(1, 1) * 4;
	my $Num = 0;
	my $Den = 0;
	
	# 分子
	$Num = GAMMAFUNCTION(($DegreesOfFreedom + 1) / 2);
	# 分母
	$Den = sqrt($DegreesOfFreedom * $Pi) * GAMMAFUNCTION($DegreesOfFreedom / 2);
	# スチューデントのt分布 Frequency Function
	$FrequencyFunction = ($Num / $Den) * ((1 + (($X * $X) / $DegreesOfFreedom)) ** (-(($DegreesOfFreedom + 1) / 2)));
	
	return $FrequencyFunction;
}
	
# 正則化された不完全ベータ関数 Regularized Incomplete Beta Function
# 引数 変数 変数 変数 ($AA, $BB, $X)
# 戻り値 正則化された不完全ベータ関数 ($RegularizedIncompleteBetaFunction)
sub REGULARIZEDINCOMPLETEBETAFUNCTION{
	my ($AA, $BB, $X) = @_;
	my $RegularizedIncompleteBetaFunction = 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;
	
	# 値の確認
	if(($AA <= 0) || ($BB <= 0) || ($X < 0) || (1 < $X)){
		return 0;
	}
	
	if($X == 0){
		# 正則化された不完全ベータ関数 Regularized Incomplete Beta Function
		$RegularizedIncompleteBetaFunction = 0;
	
		return $RegularizedIncompleteBetaFunction;
	}elsif($X == 1){
		# 正則化された不完全ベータ関数 Regularized Incomplete Beta Function
		$RegularizedIncompleteBetaFunction = 1;
	
		return $RegularizedIncompleteBetaFunction;
	}
	
	# シンプソンの公式 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;
	}
	
	# 正則化された不完全ベータ関数 Regularized Incomplete Beta Function
	$RegularizedIncompleteBetaFunction = (($SumA + $Sum + $SumB) * ($h / 3)) / &BETAFUNCTION($AA, $BB);
	
	return $RegularizedIncompleteBetaFunction;
}
	
	
# ベータ関数(近似) Beta Function Approximation
# 引数 値 ($X)
# 戻り値 ベータ関数(近似) ($BetaFunction)
sub BETAFUNCTION{
	my ($X, $Y) = @_;
	my $BetaFunction = 0;
	
	# 値の確認
	if(($X <= 0) && (abs($X - int($X)) == 0) || ($Y <= 0) && (abs($Y - int($Y)) == 0)){
		return 0;
	}
	
	# ベータ関数(近似) Beta Function Approximation
	$BetaFunction = (&GAMMAFUNCTION($X) * &GAMMAFUNCTION($Y)) / &GAMMAFUNCTION($X + $Y);
	
	return $BetaFunction;
}
	
# ガンマ関数(近似) 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;
}