FC2ブログ
最初
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/LogarithmicNormalDistribution3.txt

修正1 4/26
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/LogarithmicNormalDistribution4.txt

修正1 8/17

# 対数正規分布 Logarithmic Normal Distribution
# 引数 変数 平均 標準偏差 ($X, $Mean, $StandardDeviation)
# 戻り値 対数正規分布 (@LogarithmicNormalDistribution)
sub LOGARITHMICNORMALDISTRIBUTION{
my ($X, $Mean, $StandardDeviation) = @_;
my @LogarithmicNormalDistribution = ();

# 標準偏差の確認
if(($X < 0) || ($StandardDeviation <= 0)){
return 0;
}

# 確率密度関数 Frequency Function
$LogarithmicNormalDistribution[0] = &PROBABILITYFUNCTION($X, $Mean, $StandardDeviation);
# 下側累積確率 Lower Cumulative Distribution
$LogarithmicNormalDistribution[1] = &LOWERPROBABILITY($X, $Mean, $StandardDeviation);
# 上側累積確率 Upper Cumulative Distribution
# $LogarithmicNormalDistribution[2] = &UPPERPROBABILITY($X, $Mean, $StandardDeviation);
$LogarithmicNormalDistribution[2] = 1 - $LogarithmicNormalDistribution[1];

return @LogarithmicNormalDistribution;
}

# 確率密度 Probability Function
# 引数 変数 平均 標準偏差 ($X, $Mean, $StandardDeviation)
# 戻り値 確率密度関数 ($ProbabilityFunction)
sub PROBABILITYFUNCTION{
my ($X, $Mean, $StandardDeviation) = @_;
my $ProbabilityFunction = 0;

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

my $Pi = atan2(1, 1) * 4;
my $Ln = log($X) / log(exp(1));
my $Power = ((($Ln - $Mean) / $StandardDeviation) * (($Ln - $Mean) / $StandardDeviation));

# 対数正規分布 Frequency Function
$ProbabilityFunction = (1 / (sqrt(2 * $Pi * $StandardDeviation * $StandardDeviation) * $X)) * exp(-($Power / 2));

return $ProbabilityFunction;
}

# 下側累積確率 Lower Probability
# 引数 変数 平均 標準偏差 ($X, $Mean, $StandardDeviation)
# 戻り値 下側累積確率 ($LowerProbability)
sub LOWERPROBABILITY{
my ($X, $Mean, $StandardDeviation) = @_;
my $LowerProbability = 0;

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

my $LnX = log($X) / log(exp(1));

# 下側累積確率 Lower Probability
$LowerProbability = 0.5 + (ERRORFUNCTION(($LnX - $Mean) / sqrt(2 * ($StandardDeviation * $StandardDeviation))) / 2);

return $LowerProbability;
}

# 上側累積確率 Upper Probability
# 引数 変数 平均 標準偏差 ($X, $Mean, $StandardDeviation)
# 戻り値 上側累積確率 ($UpperProbability)
sub UPPERPROBABILITY{
my ($X, $Mean, $StandardDeviation) = @_;
my $UpperProbability = 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;

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

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

$SumA = &PROBABILITYFUNCTION($X, $Mean, $StandardDeviation);
$SumB = &PROBABILITYFUNCTION((($h * $i) - $X), $Mean, $StandardDeviation);

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

# 上側累積確率 Upper Probability
$UpperProbability = $SimpsonsRule;

return $UpperProbability;
}

# 誤差関数 Error Function
# 引数 値 ($X)
# 戻り値 誤差関数 ($ErrorFunction)
sub ERRORFUNCTION{
my ($X) = @_;
my $ErrorFunction = 0;
my $RombergsMethod = 0;
my @Romberg = ();
my $Sum = 0;
my $a = 0;
my $b = $X;
my $h = 0;
my $N = 10;
my $Power4 = 0;
my $Epsilon = 1.0e-10;
my $Pi = atan2(1, 1) * 4;

if($X == 0){
# 誤差関数 Error Function
$ErrorFunction = 0;

return $ErrorFunction;
}

# ロンバーグ積分 Romberg's Method
# (0, 0)
$Romberg[0] = (($b - $a) / 2) * (&INTEGRAND($a) + &INTEGRAND($b));

for(my $i = 1; $i < $N; $i++){
$Power4 = 1;

for(my $j = $i; $j >= 0; $j--){
if($i == $j){
$Sum = 0;
$h = ($b - $a) / (2 ** $j);

for(my $k = 1; $k <= (2 ** ($i - 1)); $k++){
$Sum += &INTEGRAND($a + (((2 * $k) - 1) * $h));
}

# ($i, 0)
$Romberg[$j] = ($Romberg[$j - 1] / 2) + ($Sum * $h);
}else {
# 累乗
$Power4 = 4 * $Power4;

# リチャードソンの補外 Richardson's Extrapolation
$Romberg[$j] = (($Power4 * $Romberg[$j + 1]) - ($Romberg[$j])) / ($Power4 - 1);

# 収束判定
if(abs($Romberg[$j] - $Romberg[$j + 1]) < $Epsilon){
$Romberg[0] = $Romberg[$j];
$i = $N;
last;
}
}
}
}

# ロンバーグ積分 Romberg's Method
$RombergsMethod = $Romberg[0];

# 誤差関数 Error Function
$ErrorFunction = (2 / sqrt($Pi)) * $RombergsMethod;

return $ErrorFunction;
}

# 被積分関数 Integrand
# 引数 変数 被積分関数 ($X, \@Integrand)
# 戻り値 被積分関数 ($Integrand)
sub INTEGRAND{
my ($X) = @_;
my $Integrand = exp(-($X * $X));

return $Integrand;
}


正規分布から対数正規分布をだす
最初
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/LogarithmicNormalDistribution.txt
修正 1
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/LogarithmicNormalDistribution2.txt

参考URL
Log-normal distribution - Wikipedia, the free encyclopedia
LOGNORMDIST - Excel - Microsoft Office Online
エクセルを用いた対数正規分布の計算
対数正規分布 - 高精度計算サイト

一言
やり方が悪いのか、正規分布から対数正規分布を出した方が早く誤差が小さい

修正1
&FREQUENCYFUNCTION($LnX, $Mean, $StandardDeviation) を $X で割った

修正1 4/26
(1 / ($StandardDeviation * $X * sqrt(2 * $Pi)))を
(1 / (sqrt(2 * $Pi * $StandardDeviation * $StandardDeviation) * $X))

修正2 8/17
誤差関数から値を出す

メモ
正規分布から出すと確率密度が違うことがある
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

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

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

検索フォーム


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