パーセンタイルの数値計算
統計学におけるパーセンタイルをプログラミングで計算します。
平均からの差 | : $\times\sigma$ |
パーセンタイル | : |
珍しさ | : |
※)$6\sigma$のとき9桁目までは正しく、$3\sigma$では16桁とも合っているようです。
「Expected Fraction of Population Inside Range」を2で割って0.5を足したものがパーセンタイルです。
正規分布関数(ガウス関数)を$\left(-\infty, \mu + s\sigma\right]$で積分します。計算の都合上、区間を$\left(-\infty, \mu\right]$と$\left[\mu, \mu + s\sigma\right]$に分け、$\left[\mu, \mu + s\sigma\right]$で積分して0.5を足します。
\begin{eqnarray*}p&=&\int_{\mu}^{\mu+s\sigma}\mathrm{d}x\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\left(x-\mu\right)^{2}}{2\sigma^2}\right)\\&=&\int_{0}^{s\sigma}\mathrm{d}x'\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{x'^{2}}{2\sigma^2}\right)\\&=&\int_{0}^{s}\mathrm{d}x''\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x''^{2}}{2}\right)\\&=&\int_{0}^{s}\mathrm{d}x''\frac{1}{\sqrt{2\pi}}\sum_{n=0}^{\infty}\frac{\left(-\frac{x''^{2}}{2}\right)^{n}}{n!}\\&=&\frac{1}{\sqrt{2\pi}}\sum_{n=0}^{\infty}\frac{\left(-\frac{1}{2}\right)^{n}}{n!}\int_{0}^{s}\mathrm{d}x''x''^{2n}\\&=&\frac{1}{\sqrt{2\pi}}\sum_{n=0}^{\infty}\frac{\left(-\frac{1}{2}\right)^{n}}{n!}\frac{1}{2n+1}\left[x''^{2n+1}\right]_{0}^{s}\\&=&\frac{1}{\sqrt{2\pi}}\sum_{n=0}^{\infty}\frac{\left(-\frac{1}{2}\right)^{n}s^{2n+1}}{n!\left(2n+1\right)}\end{eqnarray*}プログラミングでは、$\sum$を第$200$項まで計算すれば十分な精度が得られます。各項を計算するときに、精度をできるだけ落とさないために、$i$について$n$回の繰り返しを用いて、初期値を$\frac{s}{2n+1}$として、繰り返しの最中に$\frac{s^{2}}{2\times i}$を掛け、繰り返しの終了後に$\left(-1\right)^{n}$を掛けます。
この$p$に$0.5$を足して、$100$倍すれば、パーセンタイルとなります。
計算例はIQにあります。
第200項までで充分であることの根拠:
想定する最大値を$10\sigma$($s=10$)とし、パーセンタイルは小数点以下14桁まで計算するとします。
(JavaScriptの性能的に$10\sigma$は計算できませんが、マージンの意味もあります。)
和の一般項は$\displaystyle\frac{\left(-\frac{1}{2}\right)^{n}s^{2n+1}}{n!\left(2n+1\right)}$です。第$n$項と第$n+1$項の和は\[\frac{\left(-\frac{1}{2}\right)^{n}s^{2n+1}}{n!\left(2n+1\right)}+\frac{\left(-\frac{1}{2}\right)^{n+1}s^{2\left(n+1\right)+1}}{\left(n+1\right)!\left(2\left(n+1\right)+1\right)}=\frac{\left(-1\right)^ns^{2n+1}}{2^nn!\left(2n+1\right)}\left(1-\frac{s^2\left(2n+1\right)}{2\left(n+1\right)\left(2n+3\right)}\right)\]です。第$n+2$項と第$n+3$項の和は\[\frac{\left(-\frac{1}{2}\right)^{n+2}s^{2\left(n+2\right)+1}}{\left(n+2\right)!\left(2\left(n+2\right)+1\right)}+\frac{\left(-\frac{1}{2}\right)^{n+3}s^{2\left(n+3\right)+1}}{\left(n+3\right)!\left(2\left(n+3\right)+1\right)}=\frac{\left(-1\right)^ns^{2n+5}}{2^{n+2}\left(n+2\right)!\left(2n+5\right)}\left(1-\frac{s^2\left(2n+5\right)}{2\left(n+3\right)\left(2n+7\right)}\right)\]です。これらの比を計算すると、\[\frac{第n+2項と第n+3項の和}{第n項と第n+1項の和}=\frac{s^4\left(2n+1\right)}{2^2\left(n+1\right)\left(n+2\right)\left(2n+5\right)}\frac{\displaystyle\left(1-\frac{s^2\left(2n+5\right)}{2\left(n+3\right)\left(2n+7\right)}\right)}{\displaystyle\left(1-\frac{s^2\left(2n+1\right)}{2\left(n+1\right)\left(2n+3\right)}\right)}\]$n$が大きいとき、\[\frac{第n+2項と第n+3項の和}{第n項と第n+1項の和}\simeq\frac{s^4}{2^2n^2}=\left(\frac{s^2}{2n}\right)^2\]$s$が大きいほど高次の項の影響は大きいですが、$s=10$とすると、\[\frac{第n+2項と第n+3項の和}{第n項と第n+1項の和}\simeq\left(\frac{50}{n}\right)^2\]また、$n$が大きくなるにつれて高次の項の影響が小さくなります。これは1や2に比べて充分大きいときの近似ですので、$n=150$とすると、\[\frac{第n項と第n+1項の和}{\sqrt{2\pi}}=1.09\cdots\times10^{-10}\]となって、すでにパーセンタイルの小数点以下8桁ほどの影響ですが、その次の項の影響を調べると、\[\frac{第n+2項と第n+3項の和}{第n項と第n+1項の和}\simeq\left(\frac{1}{3}\right)^2\]となります。$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{3}\right)^{2n}$は$0.125$ですが、$\displaystyle\sum_{n=1}^{20}\left(\frac{1}{3}\right)^{2n}$を計算すると19桁まで一致していて、高次になるほど影響は小さくなるので$n=150+20$項目はもはや計算結果に影響しません。
なので、(結果を見つつ適当にとっていたのですが)第$200$項目までとれば充分です。
計算量の節約のために、より最高次を小さくすることも可能です。
数学的には上記の理論ですが、JavaScriptの精度的に限界があります。
$6\sigma$のとき、9桁までは合っているようです。