Pythonプログラミング（ステップ７・リスト・移動平均）
このページでは、データ処理において基本中の基本である「平均操作」について考えてみる。
平均操作
$N$個のデータ点$X_i$ $(i=0,\cdots,N-1)$ があった場合、それらの代表値として平均は $$ m = \frac{1}{N} \sum_{i=0}^{N-1} X_i \tag{1} $$ で与えられる。
各データ点に「重み」$w_i$ をつけた荷重付き平均 $$ m = \sum_{i=0}^{N-1} w_i x_i \tag{2} $$ もよく用いられる。ここで、$w_0 + w_1 + \cdots + w_{n-1} = 1$ である。 そして、(1)は荷重が一定 $w_i = 1/N$ の特殊な場合、と考えることができる。
改めて言うまでもなく、平均とは「全体的に均す」操作に他ならない。
移動平均
自然現象や社会現象によって生じる時系列は複雑に変動するため、大まかな傾向（トレンド）がなかなか掴みにくい場合が多い。 例えば、株価は目まぐるしく変化して、この先、上がるのか、下がるのか、予測が難しいし、 地球全体が温暖化していると言われてはいるものの、季節等による変動幅も大きい。
そのようなデータに対してよく使われるのが移動平均である。$i$番目の時系列点の値を$X_i$とすれば、その点での移動平均値を $$ \bar{X}_i = \frac{1}{M+1} \sum_{i-M/2}^{i+M/2} X_i \tag{3} $$ で与える（ここでは簡単のため、$M$は偶数と仮定した）。 注目する点$i$を移動させながら、その周りで平均を計算することで、「局所的な代表値」の系列を得ようとするわけである。
3点（$M=2$）での移動平均の計算
情報基礎の教材「株価と移動平均線」(学内限定）では、株価を例にとって、 複雑な時系列の傾向を把握する勘弁な手法が紹介されているので、そちらも参照すること。
以下は、東北電力が公開している東北・新潟地区の電力消費量の時系列データ（データはからダウンロードできる）を読み込み、72時間ずつの時間幅で平均して表示するコードの例である。
水色が1時間毎の時系列、青色が移動平均を取った後のデータを表す。平均を行う時間幅($m$に対応）はプログラム中で、変数 WSIZEで設定されている。
横軸は、Unix時間（1970年1月1日の深夜0時からの経過秒数）である。
データ: juyo_2018_tohoku.txt
# coding: utf-8
import pandas as pd
import matplotlib.pyplot as plt
data = pd.read_csv('juyo_2018_tohoku.txt', encoding='Shift_JIS', skiprows=1)
date = data['DATE'].str.cat(data['TIME'],sep=" ")
datetime = pd.to_datetime(date,format='%Y/%m/%d %H:%M')
unixtime = datetime.map(pd.Timestamp.timestamp)
x_raw = data['実績(万kW)'].tolist()
N = len(x_raw)
WSIZE = 72
x = []
for i in range(N):
s = 0
n = 0
for j in range(i-WSIZE//2,i+WSIZE//2+1):
if j>=0 and j<N:
s = s + x_raw[j]
n = n + 1
x.append(s/n)
plt.plot(unixtime,x_raw, color=(0.0,1.0,1.0), linewidth=1.0)
plt.plot(unixtime,x, color=(0.0,0,1.0), linewidth=2.0)
plt.xlabel('UNIXTIME')
plt.ylabel('10^4 KWH')
plt.grid(True)
plt.show()
実行結果
横軸の 1e9 は $\times 10^9$ を表す。
練習：時間幅の調整
上記のコードの平均時間幅（WSIZE）、および、プロットする時間範囲を変えて、移動平均の効果を確認してみなさい。
また、この時系列にはどのような周期的な要因が反映されているか、データに基づいて考察しなさい。
音の加工
我々の最も身近でコミュニケーションに欠かせない時系列情報は「音」であろう。 音は、空気の粗密（あるいは圧力）の変化を耳やマイクロフォンで変換したもので、コンピューター上では、一定時間間隔でのそれらの値を数値の列として表現される。
データの時間間隔は、毎秒あたりのデータ点の数（サンプリングレート、あるいはフレームレート）で表せる。 例えばCDの音源のサンプリングレートは 44100 サンプル/秒なので、1時間の音楽は、$44100 \times 3600$ （実際には左右のチャンネルがあるのでさらにその二倍）のデータ点で構成されることになる。
以下のコードは、WAV形式の音データファイル(chime.wav)を読み込み、WSIZEの移動平均を取った後、
結果を別の音データファイル(soft_chime.wav)に保存するコードの例である。
移動平均の部分は、上記の電力時系列の処理で用いたものと全く同じである。
サンプル音データ: chime.wav
WSIZEを大きくすると終了までにかなり時間がかかるので注意。
# coding: utf-8
import numpy as np
from scipy.io import wavfile
import matplotlib.pyplot as plt
data = wavfile.read('chime.wav')
framerate = data[0]
x_raw = data[1]
N = len(x_raw)
WSIZE = 16
x = np.empty((N,),dtype='int16')
frames=np.empty((N,))
for i in range(N):
s = 0
n = 0
for j in range(i-WSIZE//2,i+WSIZE//2+1):
if j>=0 and j<N:
s = s + x_raw[j]
n = n + 1
x[i]=s/n
frames[i] = i/framerate
plt.plot(frames,x_raw, color=(0.0,1.0,1.0), linewidth=0.1)
plt.plot(frames,x, color=(0.0,0,1.0), linewidth=0.1)
plt.xlabel('SEC')
plt.ylabel('AMPLITUDE')
plt.grid(True)
plt.show()
wavfile.write('soft_chime.wav',framerate,x)
Google Colabを使っている場合は、引き続いて、以下の別のセルで実行すると、加工後の音が再生できる。
from IPython.display import Audio Audio(x,rate=framerate)
移動平均の幅(WSIZE)が大きいほど、加工後の音は、高音部が聞こえにくくなる一方で、低音（lower range)には変化が無いはずである。
すなわち、移動平均は、ローパスフィルタ（低域通過フィルタ）として動作するとも言うことができる。
さらに、(2)式の荷重$w_i$に相当する値を変えることによって、さらに異なる特性のフィルターを設計することもできる（その場合、$w_i$の総和が1になるという制限は不要）。
移動平均の畳み込みによる表現
上記のコードで移動平均を計算する箇所は、効率が悪いので、NumPyの畳み込み(convolution)の機能を使うのが実際的である。
$w$のサイズが$M$(偶数)の場合は $$ \sum_{k=0}^{M} w_k \, x_{i+M/2-k-1} $$
ここで、配列$w$のサイズを$M+1$ ($M$は偶数と仮定）とすると、 配列$x$と$w$の畳み込みの結果得られる配列$y$の$i$番目の要素は $$ y_i = \sum_{k=0}^{M} w_k \, x_{i+M/2-k} $$ で定義され(ここで$M/2$は整数の除算）、NumPyではこの計算を
y = np.convolve(w, x, 'same')
と表現する。この機能を使うと、上記のサンプルコードで移動平均を求めている箇所は以下のように書き直すことができる。†
N = len(x_raw) WSIZE = 16 w = np.ones((WSIZE+1,))/(WSIZE+1) y = np.convolve(w, x_raw,'same') x = y.astype(np.int16) frames=np.linspace(0,N/framerate,N)
最後の４行のコードの流れは以下の通りである：
- まず、要素数
WSIZE+1の配列wの値を $1/(WSIZE+1)$ に設定しておく。 - 音データの入った配列
x_rawとwの畳み込みを計算し、yにセットする。 - 得られた結果を、16ビットの整数型の配列に変換し
xにセットする。 - 配列
framesは、0から $N/framerate$までの値を$N$等分した値（データ点の経過秒数）を設定する。
† 厳密には、上述の移動平均の計算コードでは、データの端点では移動平均の幅の範囲に入っている点のみの平均を計算しているのに対して、 NumPyの畳み込みでは、データ点の前後に0が続いていると仮定して計算するので、端点付近での結果には違いがある。
画像のぼかし
画像のデータは、縦横の升目状の画素（ピクセル）の集合として扱われている。一つの画素に、三原色の強度（明るさ）を対応させる場合が多い。
このピクセルデータに移動平均の考え方を適用すると、どのような効果を生むだろうか。 時系列の場合は、一次元的な「両どなり」付近の値の平均を取ったけれども、画像の場合は「上下左右」方向に 平均する必要がある点が違うだけで、データの基本的な操作には何ら変わりはない。
サンプル画像: bulbul.png
# coding: utf-8
from skimage import io,img_as_float
image1 = io.imread('bulbul.png')
image2 = img_as_float(image1)
nx = image2.shape[0]
ny = image2.shape[1]
WSIZE=4
for x in range(nx):
for y in range(ny):
sr = sg = sb = 0
n = 0
for i in range(x-WSIZE//2,x+WSIZE//2+1):
for j in range(y-WSIZE//2,y+WSIZE//2+1):
if i>=0 and i<nx and j>=0 and j<ny:
sr = sr + image1[i,j,0]
sg = sg + image1[i,j,1]
sb = sb + image1[i,j,2]
n = n + 1
image2[x,y,0] = sr/n/256
image2[x,y,1] = sg/n/256
image2[x,y,2] = sb/n/256
io.imshow(image2)
io.show()
io.imsave('bluered_bulbul.png', image2)
画像は処理すべきデータ点数が大きいので、移動平均の計算には時間を要してしまう。 音データの場合と同様に、畳込みを使って同様の処理を行うコードの例を以下に示す。 ここではSciPyライブラリの中の２次元畳み込み関数 convolve2d( )を用いている。
# coding: utf-8
from skimage import io,img_as_float
from scipy.signal import convolve2d
import numpy as np
image1 = io.imread('bulbul.png')
image2 = img_as_float(image1)
WSIZ=32
w = np.ones((WSIZ+1, WSIZ+1))/((WSIZ+1)**2)
red = convolve2d(image1[:,:,0],w,mode='same', boundary='symm')
green = convolve2d(image1[:,:,1],w,mode='same', boundary='symm')
blue = convolve2d(image1[:,:,2],w,mode='same', boundary='symm')
image2[:,:,0] = red[:,:]/256
image2[:,:,1] = green[:,:]/256
image2[:,:,2] = blue[:,:]/256
io.imshow(image2)
io.show()
io.imsave('bluer_bulbul.png', image2)
練習：画像フィルタ
直前の畳み込みを使ったサンプルプログラムで、平均の荷重を表すwを、$3 \times 3$の配列
[ [0, 1, 0], [1,-4, 1], [0, 1, 0] ]
に変えた場合、どのような画像が生成されるか試してみなさい。
ヒント
$(x,y)$の位置の画素値を$\phi(x,y)$とすると、上記の荷重に対応する演算 $$ \phi(x+1,y) + \phi(x-1,y) + \phi(x,y+1) + \phi(x,y-1) - \phi(x,y) $$ は、$\phi$の各点にラプラス演算子を施す操作に相当する。 さらに詳しくは こちらのページも参照のこと。