ソルバーで正規分布の90%範囲を求める

[PR]当サイトはアフィリエイト広告を利用しています。

はじめに
台形公式の漸増計算をするユーザー定義関数
ソルバーで90%区間の右側の値を変数にして最適化

はじめに

突然ですが次のような問題を考えます。

ある会社の健康診断で身長を測ったら、平均は170.5cm、標準偏差は5.4cmだった。

身長の分布を正規分布とみなすとき、90%のデータが入る区間を求めよ。

このとき平均から上に45%、下に45%の区間を求めるとします。

やりかたは正規分布表を用いれば簡単です。

正規分布表の値がp(u) = 0.45に最も近い位置を求めます。この場合u = 1.64です。

そしてu = (X – m) / σ

に平均mと標準偏差σとuを代入して、Xについて解けば、X = 179.36cmが求まります。

これが平均から上に45%の位置。

下に45%の位置はm – (179.36 – m) = 161.64cmと求まりました。

身長の90%は161.64cmから179.36cmの範囲にあるということです。

今回は正規分布表を用いずに、数値積分とソルバーでこの90%区間を求めてみようというお話です。

数値積分に関して参考にしたのは以下の書籍。

栗原正仁（2011）『わかりやすい数値計算入門』ムイスリ出版.

今回は台形公式の漸増計算を用いて自動的に収束判定をして分割数を自動的に決め、数値積分を用いて正規分布曲線から確率を求めます。

そしてその計算をソルバーに投げて、90%区間の右側の値を変数にして最適化します。

台形公式の漸増計算をするユーザー定義関数

台形公式の漸増計算は、積分区間[a,b]、分割数n、刻み幅h=(b-a)/nとすると

で計算できます。

今回のf(x)は

となります。f(x)は正規分布曲線を表します。

ε= 0.00001などとして

となったら収束したとみなして計算終了とします。

以上を踏まえて作成したユーザー定義関数が以下。

Function CalcSTD(ave As Variant, std As Variant, x As Variant)
Dim n As Double
Dim h As Double
Dim S1 As Double
Dim S2 As Double
Dim temp As Double
Dim e As Double
Dim i As Long
e = 0.00001
n = 1#
h = (x - ave) / n
S1 = (f(ave, ave, std) + f(x, ave, std)) * ((x - ave) / 2)
n = n + 1#
h = (x - ave) / n
S2 = S1 / 2 + h * f(ave + h, ave, std)
While (Abs(S1 - S2) > e)
n = n * 2#
h = (x - ave) / n
temp = S2
S2 = 0
For i = 1 To (CLng(n) - 1) Step 2
    S2 = S2 + f(ave + CDbl(i) * h, ave, std)
Next i
S2 = S2 * h + temp / 2
S1 = temp
Wend

CalcSTD = S2
End Function

Function f(x As Variant, ave As Variant, std As Variant)
f = (2.71828 ^ (-(x - ave) * (x - ave) / (2 * std * std))) / (((2 * 3.14159) ^ 0.5) * std)
End Function

まず31行目のFunction f()から解説します。

これは正規分布曲線のある位置での値を出力する関数です。

次に11行目のhの式。

積分区間は[m, X]なのでこのコードになります。

12行目から15行目は最初のSnとSn/2を求めています。

もちろん収束していないので16行目から26行目のWhile文に続きます。

21行目から23行目のFor文でのCLng(n)とCDbl(i)はiとnをLongとDoubleに型変換しています。

19行目から24行目でSnを求めて、S1にSn/2を格納しています。

このWhile文で収束するまで分割数nを増やしていき、収束したら終了します。

28行目でS2をユーザー定義関数の出力とします。