標準正規分布の2.5%上側分位点が約1.96であることはどうやって計算するのかという質問を受けた． 分位点は累積分布関数の逆関数なので，まず累積分布関数 を計算し，その逆関数を求める計算になる (世の中のソフトウェアはもっと能率的な実装が取られるはずだが，それはひとまず置く)．

信頼区間の計算などでは分位点の検索では裾(上側2.5%,1%,0.5%など)を調べることが多いので，まずは 裾の確率の計算を考える． \[ Q(x) = \int_x^\infty \frac{1}{2\pi} \exp\left(-x^2/2\right) dx \] は確率変数\(X\sim N(0,\sigma^2)\)が\(X > x\)を取る確率であり，例えば\(Q(1.96) \approx 2.5\%\)である ことが知られている．これを自力で計算したい．

\( t = x^2/2 \) とおくと， \[ Q(x) = \frac{1}{4\pi} \int_{x^2/2}^\infty \frac{ e^{-t} }{\sqrt{t} } dt \] なので，この問題は積分 \[ f(x) = \int_x^\infty t^{-1/2} e^{-t} dt \] の計算に帰着される．これは\(t^{-1/2}\)という\(t\to\infty\)で 0に収束する因子が掛かっているが，部分積分を繰り返すことで，さらに早い収束をする因子に書き換えられる．

\begin{eqnarray*} f(x) &=& \int_x^\infty t^{-1/2} e^{-t} dt \\ & = & \left. -e^{-t}t^{-1/2} \right|_x^\infty - \frac{1}{2} \int_x^\infty t^{-3/2} e^{-t} dt \\ & = & e^{-x}x^{-1/2} - \frac{1}{2} \int_x^\infty t^{-3/2} e^{-t} dt \\ & = & e^{-x}x^{-1/2} + \left. \frac{1}{2} e^{-t}t^{-3/2} \right|_x^\infty + \frac{1\cdot 3}{2\cdot 2} \int_x^\infty t^{-5/2} e^{-t} dt \\ & = & e^{-x}x^{-1/2} - \frac{1}{2} e^{-x}x^{-3/2} + \frac{1\cdot 3}{2\cdot 2} \int_x^\infty t^{-5/2} e^{-t} dt \\ & = & e^{-x}x^{-1/2} - \frac{1}{2} e^{-x}x^{-3/2} + e^{-x}x^{-1/2} - \left. \frac{1 \cdot 3}{2 \cdot 2} e^{-t}t^{-5/2} \right|_x^\infty - \frac{1\cdot 3 \cdot 5}{2 \cdot 2 \cdot 2 } \int_x^\infty t^{-7/2} e^{-t} dt \end{eqnarray*} この操作は無限に続けられる．整理すると \begin{eqnarray*} f(x) &=& \frac{e^{-x}}{\sqrt{x}} \left\{ 1 - \frac{1}{2x} + \frac{3}{4x^2} - \frac{15}{8 x^3} + \frac{105}{16x^4} - \cdots \right\} \\ Q(x) &=& f(x^2/2)/\sqrt{4\pi} \end{eqnarray*}

中括弧の中を3次で打ち切った近似とRによる高精度の値との比較する表を以下に示す．結果を見ると，この展開は\(x > 3\)で有効と判断できる．なお，この級数は発散するもので，\(xの値に応じて\)高次の項まで取りすぎると精度がかえって悪化する．

この計算には以下のRスクリプトを使った．

f4 = function(x) exp(-x)/sqrt(x)*( 1 - 1/(2*x) + 3/(4*x*x));
Q4 = function(a) f4(a*a/2)/sqrt(4*pi) ;
xs = c(1,1.96,2,3,5,10)
tab = cbind(xs,Q4(xs),pnorm(xs,lower.tail=FALSE))