FC2ブログ

# ベータ分布 (逆関数) Beta Distribution (Inverse Function)
# 引数 # 引数 変数 変数 変数 ($X, $A, $B)
# 戻り値 ベータ分布 (逆関数) (@InverseFunction)
sub BETADISTRIBUTIONINVERSE{
my ($X, $A, $B) = @_;
my @InverseFunction = ();

# 累積確率 変数の確認
if(($X < 0) || (1 < $X) || ($A <= 0) || ($B <= 0)){
return "Error";
}

# ベータ分布 (逆関数) Beta Distribution (Inverse Function)
# 下側 Lower
$InverseFunction[0] = &BISECTIONMETHOD($X, $A, $B);
# 上側 Upper
$InverseFunction[1] = &BISECTIONMETHOD((1 - $X), $A, $B);

return @InverseFunction;
}

# 二分法 Bisection Method
# 引数 累積確率 変数 変数 ($X, $A, $B)
# 戻り値 二分法 ($BisectionMethod)
sub BISECTIONMETHOD{
my ($X, $A, $B) = @_;
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((&BETADISTRIBUTION($X1, $A, $B) - $X) > 0);
$X2++ while((&BETADISTRIBUTION($X2, $A, $B) - $X) < 0);

# 初期値
$F_x1 = &BETADISTRIBUTION($X1, $A, $B) - $X;
$F_x2 = &BETADISTRIBUTION($X2, $A, $B) - $X;

# 計算
for(my $i = 0; $i < $Limit; $i++){
# 一つ前
$PrevMiddle = $Middle;
# 中間点
$Middle = ($X1 + $X2) / 2;

# f(Middle)
$F_m = &BETADISTRIBUTION($Middle, $A, $B) - $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;
}

# ベータ分布 Beta Distribution
# 引数 変数 変数 変数 ($X, $A, $B)
# 戻り値 ベータ分布 ($BetaDistribution)
sub BETADISTRIBUTION{
my ($X, $A, $B) = @_;
my $BetaDistribution = 0;

# 下側累積確率 Lower Probability
$BetaDistribution = &REGULARIZEDINCOMPLETEBETAFUNCTION($A, $B, $X);

return $BetaDistribution;
}

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

# 変数の確認
if(($X < 0) || (1 < $X)){
return 0;
}

if(($X == 0) || ($X == 1)){
# 正則化された不完全ベータ関数 Regularized Incomplete Beta Function
$RegularizedIncompleteBetaFunction = ($X == 0 ? 0 : 1);
}else {
# 正則化された不完全ベータ関数 Regularized Incomplete Beta Function
$RegularizedIncompleteBetaFunction = &INCOMPLETEBETAFUNCTION($AA, $BB, $X) / &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;
}

# 不完全ベータ関数 Incomplete Beta Function
# 引数 変数 変数 変数 ($AA, $BB, $X)
# 戻り値 不完全ベータ関数 ($IncompleteBetaFunction)
sub INCOMPLETEBETAFUNCTION{
my ($AA, $BB, $X) = @_;
my $IncompleteBetaFunction = 0;
my $TrapezoidalRule = 0;
my $T = 0;
my $Sum = 0;
my $HyperbolicSine = 0;
my $HyperbolicCosine = 0;
my $TempA = 0;
my $TempB = 0;
my $a = 0;
my $b = $X;
my $N = 100;
my $Pi = atan2(1, 1) * 4;
my $Pi_2 = $Pi / 2;
my $BA_2 = (($b - $a) / 2);
my $h = (2 * $Pi) / $N;
my $Begin = -($N / 2);
my $End = $N / 2;

# 値の確認
if(($AA <= 0) || ($BB <= 0) || ($X < 0) || (1 < $X)){
return 0;
}

# 計算
for(my $i = $Begin; $i <= $End; $i++){
$HyperbolicSine = &HYPERBOLICSINE($i * $h);
$HyperbolicCosine = &HYPERBOLICCOSINE($Pi_2 * $HyperbolicSine);
$T = (($BA_2 * &HYPERBOLICTANGENT($Pi_2 * $HyperbolicSine)) + (($b + $a) / 2));
$TempA = &FUNCTION($AA, $BB, $T);
$TempB = &HYPERBOLICCOSINE($i * $h) / ($HyperbolicCosine * $HyperbolicCosine);

$Sum += ($TempA * $TempB);
}

# 台形公式 Trapezoidal Rule
$TrapezoidalRule = ((($b - $a) * $Pi * $h) / 4) * $Sum;

# 不完全ベータ関数 Incomplete Beta Function
$IncompleteBetaFunction = $TrapezoidalRule;

return $IncompleteBetaFunction;
}

# 関数 Function
# 引数 積分変数 被積分関数 (\@Integrand)
# 戻り値 関数 ($Function)
sub FUNCTION{
my ($AA, $BB, $X) = @_;
my $Function = ($X ** ($AA - 1)) * ((1 - $X) ** ($BB - 1));

return $Function;
}

# 双曲線正弦 Hyperbolic Sine
# 引数 値 ($X)
# 戻り値 双曲線正弦 ($HyperbolicSine)
sub HYPERBOLICSINE{
my ($X) = @_;
my $HyperbolicSine = 0;

# 双曲線正弦 Hyperbolic Sine
$HyperbolicSine = (exp($X) - exp(-$X)) / 2;

return $HyperbolicSine;
}

# 双曲線余弦 Hyperbolic Cosine
# 引数 値 ($X)
# 戻り値 双曲線余弦 ($HyperbolicCosine)
sub HYPERBOLICCOSINE{
my ($X) = @_;
my $HyperbolicCosine = 0;

# 双曲線余弦 Hyperbolic Cosine
$HyperbolicCosine = (exp($X) + exp(-$X)) / 2;

return $HyperbolicCosine;
}

# 双曲線正接 Hyperbolic Tangent
# 引数 値 ($X)
# 戻り値 双曲線正接 ($HyperbolicTangent)
sub HYPERBOLICTANGENT{
my ($X) = @_;
my $HyperbolicSine = (exp($X) - exp(-$X)) / 2;
my $HyperbolicCosine = (exp($X) + exp(-$X)) / 2;
my $HyperbolicTangent = 0;

# 双曲線正接 Hyperbolic Tangent
$HyperbolicTangent = $HyperbolicSine / $HyperbolicCosine;

return $HyperbolicTangent;
}


参考URL
ベータ分布(逆関数) - 高精度計算サイト
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

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

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

検索フォーム


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

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。