# スチューデントのt分布 (逆関数) Student's T-Distribution (Inverse Function)
# 引数 累積確率 自由度 ($X, $DegreesOfFreedom)
# 戻り値 スチューデントのt分布 (逆関数) (@InverseFunction)
sub STUDENTSTDISTRIBUTIONINVERSE{
	my ($X, $DegreesOfFreedom) = @_;
	my @InverseFunction = ();
	
	# 累積確率 標準偏差の確認
	if(($X < 0) || (1 < $X) || ($DegreesOfFreedom <= 0)){
		return "Error";
	}
	
	if(($X == 0) || ($X == 1)){
		return "Inf or 0";
	}
	
	# スチューデントのt分布 (逆関数) Student's T-Distribution (Inverse Function)
	# 下側 Lower
	$InverseFunction[0] = &BISECTIONMETHOD($X, $DegreesOfFreedom);
	# 上側 Upper
	$InverseFunction[1] = &BISECTIONMETHOD((1 - $X), $DegreesOfFreedom);
	
	return @InverseFunction;
}
	
# 二分法 Bisection Method
# 引数 累積確率 平均 標準偏差 ($X, $DegreesOfFreedom)
# 戻り値 二分法 ($BisectionMethod)
sub BISECTIONMETHOD{
	my ($X, $DegreesOfFreedom) = @_;
	my $BisectionMethod = 0;
	my $X1 = 0;
	my $X2 = 0;
	my $F_m = 0;
	my $F_x1 = 0;
	my $F_x2 = 0;
	my $Middle = 0;
	my $PrevMiddle = 0;
	my $Limit = 100;
	my $Epsilon = 1.0e-20;
	
	# 区間
	$X1-- while((&STUDENTSTDISTRIBUTION($X1, $DegreesOfFreedom) - $X) > 0);
	$X2++ while((&STUDENTSTDISTRIBUTION($X2, $DegreesOfFreedom) - $X) < 0);
	
	# 初期値
	$F_x1 = &STUDENTSTDISTRIBUTION($X1, $DegreesOfFreedom) - $X;
	$F_x2 = &STUDENTSTDISTRIBUTION($X2, $DegreesOfFreedom) - $X;
	
	# 計算
	for(my $i = 0; $i < $Limit; $i++){
		# 一つ前
		$PrevMiddle = $Middle;
		# 中間点
		$Middle = ($X1 + $X2) / 2;
	
		# f(Middle)
		$F_m = &STUDENTSTDISTRIBUTION($Middle, $DegreesOfFreedom) - $X;
	
		# 置き換え
		if(($F_m * $F_x1) > 0){
			$X1 = $Middle
		}
		elsif(($F_m * $F_x2) > 0){
			$X2 = $Middle;
		}
	
		# 二分法 Bisection Method
		$BisectionMethod = $Middle;
	
		# 収束判定
		last if(abs($Middle - $PrevMiddle) < $Epsilon);
	}
	
	return $BisectionMethod;
}
	
# スチューデントのt分布 Student's T-Distribution
# 引数 変数 自由度 ($X, $DegreesOfFreedom)
# 戻り値 スチューデントのt分布 ($StudentsTDistribution)
sub STUDENTSTDISTRIBUTION{
	my ($X, $DegreesOfFreedom) = @_;
	my $StudentsTDistribution = 0;
	my $t = ($X + sqrt(($X * $X) + $DegreesOfFreedom)) / (2 * sqrt(($X * $X) + $DegreesOfFreedom));
	
	# 下側累積確率 Lower Probability
	$StudentsTDistribution = &REGULARIZEDINCOMPLETEBETAFUNCTION(($DegreesOfFreedom / 2), ($DegreesOfFreedom / 2), $t);
	
	return $StudentsTDistribution;
}
	
# 正則化された不完全ベータ関数 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;
}