一寸先は闇
FC2ブログ

一寸先は闇

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

最新記事

perl で特殊関数 40 対数積分 li(x) (Logarithmic Integral li: LOGARITHMICINTEGRALLI)

2010 09/11


# 対数積分 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) - 高精度計算サイト
Category : 

perl

 > 

perlで特殊関数

Tags : 

perl 特殊関数 対数積分li(x) LOGARITHMICINTEGRALLI 

perl で特殊関数 39 指数積分 Ei(x) (Exponential Integral Ei: EXPONENTIALINTEGRALEI)

2010 09/04


# 指数積分 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) - 高精度計算サイト
Category : 

perl

 > 

perlで特殊関数

Tags : 

perl 特殊関数 指数積分 Ei(x) EXPONENTIALINTEGRALEI 

perl で特殊関数 38 指数積分 En(x) (Exponential Integral En: EXPONENTIALINTEGRALEN)

2010 09/03


# 指数積分 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の場合は、このやり方でいいのかな？
Category : 

perl

 > 

perlで特殊関数

Tags : 

perl 特殊関数 指数積分 En(x) EXPONENTIALINTEGRALEN 

perlで統計 34 ロジスティック分布 (逆関数) (Logistic Distribution (Inverse Function): LOGISTICDISTRIBUTIONINVERSE)

2010 09/02


# ロジスティック分布 (逆関数) 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
ロジスティック分布（逆関数） - 高精度計算サイト
Category : 

perl

 > 

perlで統計

Tags : 

perl 統計 ロジスティック分布(逆関数) LOGISTICDISTRIBUTIONINVERSE 

perlで統計 33 コーシー分布 (逆関数) (Cauchy Distribution (Inverse Function): CAUCHYDISTRIBUTIONINVERSE)

2010 09/02


# コーシー分布 (逆関数) 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
コーシー分布（逆関数） - 高精度計算サイト
Category : 

perl

 > 

perlで統計

Tags : 

perl 統計 コーシー分布(逆関数) CAUCHYDISTRIBUTIONINVERSE 

次のページへ
シリーズ
月別アーカイブ カテゴリ
オンライン コンパイラ/インタプリタ
Perl 入門
Perl リファレンス
PHP 入門
PHP リファレンス
PHP
株式
テクニカル分析 計算式
テクニカル分析
数学
計算式ライブラリ
プロフィール

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

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

検索フォーム


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