一寸先は闇

手探りながらも、ごり押しでPerlのプログラムを書いてます。

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

# 対数積分 li(x) Logarithmic Integral li
# 引数 変数 ($X)
# 戻り値 対数積分 En(x) ($LogarithmicIntegralli)
sub LOGARITHMICINTEGRALLI{
my ($X) = @_;
my $LogarithmicIntegralli = 0;

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

if($X == 0){
# 対数積分 li(x) Logarithmic Integral li
$LogarithmicIntegralli = 0;
}
elsif($X == 1){
# 対数積分 li(x) Logarithmic Integral li
$LogarithmicIntegralli = "Inf";
}else {
# 対数積分 li(x) Logarithmic Integral li
$LogarithmicIntegralli = &EXPONENTIALINTEGRALEI(log($X) / log(exp(1)));
}

return $LogarithmicIntegralli;
}

# 指数積分 Ei(x) Exponential Integral Ei
# 引数 変数 ($X)
# 戻り値 指数積分 En(x) ($ExponentialIntegralEi)
sub EXPONENTIALINTEGRALEI{
my ($X) = @_;
my $ExponentialIntegralEi = 0;
my $EulersConstant = 0.577215665;
my $SimpsonsRule = 0;
my $h = 0.1;
my $i = 1;
my $Sum = 0;
my $SumA = 0;
my $SumB = 0;
my $PrevSum = 0;
my $Factorial = 1;
my $Limit = 100;
my $Epsilon = 1.0e-20;

if($X == 0){
# 指数積分 En(x) Exponential Integral
$ExponentialIntegralEi = "-Inf";

return $ExponentialIntegralEi;
}

if(0 < $X){
for(my $i = 1; $i < $Limit; $i++){
# 階乗
$Factorial *= $i;

$Sum += (($X ** $i) / ($i * $Factorial));
}

# 指数積分 Ei(x) Exponential Integral Ei
$ExponentialIntegralEi = $EulersConstant + (log($X) / log(exp(1))) + $Sum;
}else {
do {
# 刻み幅 $h は適当
my $T = -$X + ($i * $h);
my $Number = ($i % 2 == 0 ? 2 : 4);
my $Temp = &INTEGRAND($T);

$PrevSum = $Sum;
$Sum += $Number * $Temp;
$i++;
} while(($Sum - $PrevSum) >= $Epsilon);

$SumA = &INTEGRAND(-$X);
$SumB = &INTEGRAND((($h * $i) - $X));

# シンプソンの公式 Simpson's Rule
$SimpsonsRule = ($SumA + $Sum + $SumB) * ($h / 3);

# 指数積分 Ei(x) Exponential Integral Ei
$ExponentialIntegralEi = -$SimpsonsRule;
}

return $ExponentialIntegralEi;
}

# 被積分関数 Integrand
# 引数 変数 ($T)
# 戻り値 被積分関数 ($Integrand)
sub INTEGRAND{
my ($T) = @_;
my $Integrand = 0;

# 被積分関数 Integrand
$Integrand = exp(-$T) / $T;

return $Integrand;
}


参考URL
対数積分 - Wikipedia
対数積分 li(x) - 高精度計算サイト

# 指数積分 Ei(x) Exponential Integral Ei
# 引数 変数 ($X)
# 戻り値 指数積分 En(x) ($ExponentialIntegralEi)
sub EXPONENTIALINTEGRALEI{
my ($X) = @_;
my $ExponentialIntegralEi = 0;
my $EulersConstant = 0.577215665;
my $SimpsonsRule = 0;
my $h = 0.1;
my $i = 1;
my $Sum = 0;
my $SumA = 0;
my $SumB = 0;
my $PrevSum = 0;
my $Factorial = 1;
my $Limit = 100;
my $Epsilon = 1.0e-20;

if($X == 0){
# 指数積分 En(x) Exponential Integral
$ExponentialIntegralEi = "-Inf";

return $ExponentialIntegralEi;
}

if(0 < $X){
for(my $i = 1; $i < $Limit; $i++){
# 階乗
$Factorial *= $i;

$Sum += (($X ** $i) / ($i * $Factorial));
}

# 指数積分 Ei(x) Exponential Integral Ei
$ExponentialIntegralEi = $EulersConstant + (log($X) / log(exp(1))) + $Sum;
}else {
do {
# 刻み幅 $h は適当
my $T = -$X + ($i * $h);
my $Number = ($i % 2 == 0 ? 2 : 4);
my $Temp = &INTEGRAND($T);

$PrevSum = $Sum;
$Sum += $Number * $Temp;
$i++;
} while(($Sum - $PrevSum) >= $Epsilon);

$SumA = &INTEGRAND(-$X);
$SumB = &INTEGRAND((($h * $i) - $X));

# シンプソンの公式 Simpson's Rule
$SimpsonsRule = ($SumA + $Sum + $SumB) * ($h / 3);

# 指数積分 Ei(x) Exponential Integral Ei
$ExponentialIntegralEi = -$SimpsonsRule;
}

return $ExponentialIntegralEi;
}

# 被積分関数 Integrand
# 引数 変数 ($T)
# 戻り値 被積分関数 ($Integrand)
sub INTEGRAND{
my ($T) = @_;
my $Integrand = 0;

# 被積分関数 Integrand
$Integrand = exp(-$T) / $T;

return $Integrand;
}


参考URL
指数積分 - Wikipedia
指数積分Ei(x) - 高精度計算サイト

# 指数積分 En(x) Exponential Integral En
# 引数 変数 変数 ($N, $X)
# 戻り値 指数積分 En(x) ($ExponentialIntegralEn)
sub EXPONENTIALINTEGRALEN{
my ($N, $X) = @_;
my $ExponentialIntegralEn = 0;
my $SimpsonsRule = 0;
my $h = 0.1;
my $i = 1;
my $Sum = 0;
my $SumA = 0;
my $SumB = 0;
my $PrevSum = 0;
my $Epsilon = 1.0e-20;

if($X == 0){
my $Temp = 0;

if($N <= 1){
$Temp = "Inf";
}else {
$Temp = (($N - int($N)) == 0 ? (1 / ($N - 1)) : 0);
}

# 指数積分 En(x) Exponential Integral
$ExponentialIntegralEn = $Temp;

return $ExponentialIntegralEn;
}

do {
# 刻み幅 $h は適当
my $T = 1 + ($i * $h);
my $Number = ($i % 2 == 0 ? 2 : 4);
my $Temp = &INTEGRAND($N, $X, $T);

$PrevSum = $Sum;
$Sum += $Number * $Temp;
$i++;
} while(($Sum - $PrevSum) >= $Epsilon);

$SumA = &INTEGRAND($N, $X, 1);
$SumB = &INTEGRAND($N, $X, (($h * $i) - $X));

# シンプソンの公式 Simpson's Rule
$SimpsonsRule = ($SumA + $Sum + $SumB) * ($h / 3);

# 指数積分 En(x) Exponential Integral
$ExponentialIntegralEn = $SimpsonsRule;

return $ExponentialIntegralEn;
}

# 被積分関数 Integrand
# 引数 変数 変数 変数 ($N, $X, $T)
# 戻り値 被積分関数 ($Integrand)
sub INTEGRAND{
my ($N, $X, $T) = @_;
my $Integrand = 0;

# 被積分関数 Integrand
$Integrand = exp(-($X * $T)) / ($T ** $N);

return $Integrand;
}


参考URL
指数積分 - Wikipedia
指数積分En(x) - 高精度計算サイト

一言
Xが0の場合は、このやり方でいいのかな?

# ロジスティック分布 (逆関数) Logistic Distribution (Inverse Function)
# 引数 # 引数 変数 変数 変数 ($X, $A, $B)
# 戻り値 ロジスティック分布 (逆関数) (@InverseFunction)
sub LOGISTICDISTRIBUTIONINVERSE{
my ($X, $A, $B) = @_;
my @InverseFunction = ();

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

if(($X == 0) || ($X == 1)){
return "-Inf or Inf";
}

# ロジスティック分布 (逆関数) Logistic 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((&LOGISTICDISTRIBUTION($X1, $A, $B) - $X) > 0);
$X2++ while((&LOGISTICDISTRIBUTION($X2, $A, $B) - $X) < 0);

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

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

# f(Middle)
$F_m = &LOGISTICDISTRIBUTION($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;
}

# ロジスティック分布 Logistic Distribution
# 引数 変数 変数 変数 ($X, $A, $B)
# 戻り値 ロジスティック分布 ($LogisticDistribution)
sub LOGISTICDISTRIBUTION{
my ($X, $A, $B) = @_;
my $LogisticDistribution = 0;

# 下側累積確率 Lower Probability
$LogisticDistribution = 1 / (1 + exp(-(($X - $A) / $B)));

return $LogisticDistribution;
}


参考URL
ロジスティック分布(逆関数) - 高精度計算サイト

# コーシー分布 (逆関数) Cauchy Distribution (Inverse Function)
# 引数 # 引数 変数 変数 変数 ($X, $A, $B)
# 戻り値 コーシー分布 (逆関数) (@InverseFunction)
sub CAUCHYDISTRIBUTIONINVERSE{
my ($X, $A, $B) = @_;
my @InverseFunction = ();

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

if(($X == 0) || ($X == 1)){
return "-Inf or Inf";
}

# コーシー分布 (逆関数) Cauchy 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((&CAUCHYDISTRIBUTION($X1, $A, $B) - $X) > 0);
$X2++ while((&CAUCHYDISTRIBUTION($X2, $A, $B) - $X) < 0);

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

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

# f(Middle)
$F_m = &CAUCHYDISTRIBUTION($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;
}

# コーシー分布 Cauchy Distribution
# 引数 変数 変数 変数 ($X, $A, $B)
# 戻り値 コーシー分布 ($CauchyDistribution)
sub CAUCHYDISTRIBUTION{
my ($X, $A, $B) = @_;
my $CauchyDistribution = 0;
my $Pi = atan2(1, 1) * 4;

# 下側累積確率 Lower Probability
$CauchyDistribution = (atan2(($X - $A), $B) / $Pi) + 0.5;

return $CauchyDistribution;
}


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

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

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

検索フォーム


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

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