はじめに
Greville（グレヴィル）の平滑化とは、系列データの平滑化アルゴリズムの一種である。例えば、年齢別の死亡率$q_x$は、本来は何らかの規則的な法則に従っていると考えられる。しかしながら、実際にデータを用いて計測すると、サンプルサイズが有限であることなどの理由で、推計値には偶然的変動が含まれることになる。その偶然的変動をできる限り取り除く手段として、平滑化は有用である。実際、Grevilleの平滑化は、（簡易・完全）生命表や生保標準生命表の作成に用いられている。
・生命表の概念・作成方法 （厚労省）
・標準生命表２０１８の作成過程(アクチュアリー会)
Grevilleの平滑化の手続き
系列データ$y_{0},y_1,\cdots,y_{N-1}$が与えられているとする。まず、外挿値$y_{-1},\cdots,y_{-m}$と$y_{N},\cdots,y_{N+m}$を、それぞれ以下の式で計算する。
\begin{align}
y_{i} &= \sum_{t=1}^m a_t y_{i+t}\\
y_{i} &= \sum_{t=1}^m a_t y_{i-t}
\end{align}
そして、$i=0,1,\cdots,m-1$に対して$\tilde y_i$を以下の式で計算する
\tilde y _i = \sum_{t=-m}^m L_ty_{i+t}
このようにして得られた$\tilde y_i$が平滑化されたデータである。
上記で用いた係数$a_t$と$L_t$はいくつかのタイプがあり、それぞれ「$r$次$(2m+1)$公式」と呼ばれる。この他にも$k$で表されるパラメータがある（下記「係数の特徴付け」を参照。）。主に$r=3,k=3$が使われるので、それらの係数を以下に掲載しておく。
[Grevilleの平滑化の係数表]
◎ 3次5項式 ($k=3,r=3,m=2$)
| t | 0 | 1 | 2 |
|---|---|---|---|
| $L_{t}$ | 0.559440 | 0.293706 | -0.073426 |
| $a_t$ | - | 2.000000 | -1.000000 |
◎ 3次7項式 ($k=3,r=3,m=3$)
| t | 0 | 1 | 2 | 3 |
|---|---|---|---|---|
| $L_t$ | 0.412588 | 0.293706 | 0.058741 | -0.058741 |
| $a_t$ | - | 1.618034 | -0.236068 | -0.381966 |
◎ 3次9項式 ($k=3,r=3,m=4$)
| t | 0 | 1 | 2 | 3 | 4 |
|---|---|---|---|---|---|
| $L_t$ | 0.331140 | 0.266557 | 0.118470 | -0.009873 | -0.040724 |
| $a_t$ | - | 1.352613 | 0.114696 | -0.287231 | -0.180078 |
◎ 3次11項式 ($k=3,r=3,m=5$)
| t | 0 | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|---|
| $L_t$ | 0.277944 | 0.238693 | 0.141268 | 0.035723 | -0.026792 | -0.027864 |
| $a_t$ | - | 1.160811 | 0.281079 | -0.140968 | -0.204545 | -0.096377 |
◎ 3次13項式 ($k=3,r=3,m=6$)
| t | 0 | 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|---|---|
| $L_t$ | 0.240058 | 0.214337 | 0.147356 | 0.065492 | 0.000000 | -0.027864 | -0.019350 |
| $a_t$ | - | 1.016301 | 0.360880 | -0.021625 | -0.160909 | -0.138330 | -0.056317 |
（※1）負のインデックスは$L_{-t}=L_t$である。
（※2）Greville(1979)より転記した。
サンプルコード(Python)
import numpy as np
import matplotlib.pyplot as plt
#サンプルデータ作成
N = 100
y = 1+np.sin(np.arange(N)/N)+0.05*np.sin(200*np.arange(N)/N)
#Greville補正のパラメータ設定 3次9項式 (k=3,r=3,m=4)
m = 4
a = np.array([0.000000,1.352613,0.114696,-0.287231,-0.180078])
L = np.array([
-0.040724,
-0.009873,
0.11847,
0.266557,
0.33114,
0.266557,
0.11847,
-0.009873,
-0.040724
])
#外挿値の計算
y_L = np.zeros(m)
y_R = np.zeros(m)
y_ext = np.hstack([y_L,y,y_R])
for t in range(m):
y_ext[m-t-1] = np.sum(a*y_ext[m-t-1:m-t-1+m+1])
y_ext[m+N+t] = np.sum(np.flip(a)*y_ext[m+N+t-(m+1):m+N+t])
#平滑化処理
v = np.zeros(y_ext.shape)
for t in range(m,m+N):
v[t] = np.sum(L*y_ext[t-m:t+m+1])
v = v[m:m+N]
#確認
plt.plot(y)
plt.plot(v)
係数の特徴付け
係数$a_t$と$L_t$の値はどのようにして得られたのであろうか。
まず、$L_t$については、高々$r$次多項式に従う系列を保存する（＝さらなる平滑化をしない）という制約条件のもとで、平滑後の系列の$k$階微分が平均的に最小化されるように選んでいる。
また、外挿値を計算するための$a_t$は、簡単に説明するのは難しいが、$L_t$による平滑化が「最も自然に」作用するような外挿が選んでいることになる。
これらの詳細については、次回の記事で解説することにする。
参考文献
[1979] T.N.E. Greville. Moving-weighted-average smoothing extended to the extremities of tha data, MRC Technical Summary Report # 2025. Mathematics Research Center, University of Wisconsin.