Grevilleの平滑化について


# はじめに
　Greville（グレヴィル）の平滑化とは、系列データの平滑化アルゴリズムの一種である。例えば、年齢別の死亡率$q_x$は、本来は何らかの規則的な法則に従っていると考えられる。しかしながら、実際にデータを用いて計測すると、サンプルサイズが有限であることなどの理由で、推計値には偶然的変動が含まれることになる。その偶然的変動をできる限り取り除く手段として、平滑化は有用である。実際、Grevilleの平滑化は、（簡易・完全）生命表や生保標準生命表の作成に用いられている。

・[生命表の概念・作成方法 （厚労省）](https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00450012&tstat=000001031336&cycle=7&tclass1=000001060864&tclass2=000001163166&tclass3val=0)
・[標準生命表２０１８の作成過程(アクチュアリー会)](https://www.actuaries.jp/lib/standard-life-table/pdf/seimeihyo2018-katei.pdf)

# Grevilleの平滑化の手続き

　系列データ$y_{0},y_1,\cdots,y_{N-1}$が与えられているとする。まず、外挿値$y_{-1},\cdots,y_{-m}$と$y_{N},\cdots,y_{N+m}$を、それぞれ以下の式で計算する。

```math 
\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$を以下の式で計算する

```math 
\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)

```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$による平滑化が「最も自然に」作用するような外挿が選んでいることになる。
　これらの詳細については、[次回の記事](https://qiita.com/yuska/items/e92a1653e4da26c81c7f)で解説することにする。


# 参考文献
[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.