カーネル密度推定グラフの描画 with Excel

イントロダクション

Step 0

下表はある調査により得た1変数の標本です(n=30)。

カーネル密度推定・元データ

この標本をヒストグラムに通してみると，二峰性の分布であることがわかります。

こうした特徴を持つ標本で，カーネルカーブを描画したいと思います。以下，Excelでこれをおこなうためのstepです。

作図にあたり参考にした資料

工程

Step 1前提

カーネル密度推定量f̂(x)は下式にて求めます。

$[LaTeX] \displaystyle \hat{f}_h(x) = \frac{1}{nh} \sum_{i=1}^{n} K \left( \frac{x-x_i}{h} \right)$

ここで前提として，「R」のdensity()関数デフォルトの設定，すなわち引数bwをnrd0でプロットしたときと同様の成果物を得ることを目的に据えたいと思います。すなわちバンド幅hは，Silvermanの方法

$[LaTex] \displaystyle h = \frac{0.9 \times min \Bigl\{ s, \frac{IQR}{1.34} \Bigr\} }{n^{\frac{1}{5}}}$

(ただしsは標準偏差)で求め，カーネル関数は

$[LaTex] \displaystyle K(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}$

すなわちガウス関数を利用するものとします。

以下これらの式は，順に❶式，❷式，❸式と呼称します。

Step 2「2」式の計算(1)

❷式をシートの上方で求めます。

とはいえ現状では当該部分に利用可能なスペースが残されていないので，シートの最上部にあらたに3つの行を追加したうえで，

ここ(セルA1, A2)に見出し「n」「h」を用意することにします。

見出しn,hの追加

Step 3「2」式の計算(2)

「n」については関数で求めるか，直接タイプしておきます(下式は“関数”の場合)。

つづいて「h」を求めます。これは，❷式をExcel上の計算式にて配置します。

B1	=COUNT(A5:A34)
B2	=( 0.9 * MIN( STDEV.S(A5:A34), (QUARTILE.INC(A5:A34, 3)-QUARTILE.INC(A5:A34, 1)) /1.34 ) ) /B1^(1/5)

n,hを求める

Step 4「3」式の計算

この作例でいう変数名「x」の直右の列から，❶式にいういくらかのxを用意します。換言すれば，成果物としたいグラフ横軸のレンジをイメージして，この数字を手繰っておく必要があります。

具体的に，ここでは0－100の範囲でグラフを描きたいと考えています。ということで0から100までを任意のスパンで埋めていく必要があるのですが，当然，より短いスパンのほうがより忠実なグラフを描画することができます。もっとも，その分処理すべき計算の量がかさむこととなるので，どの程度の“忠実”を目指すかで適宜斟酌すればいいかと思います。

ともあれ，この作例では下図のように0－100を増分2で埋めることに決めました。

0-100まで増分2でxを作成