numpyでconvolveしてみる

numpyを使って数値計算で畳み込みをしてみたのでメモしておきます。

numpyで畳み込みするにはnumpy.convolveという関数を使用します。
\[f(x) = x^2(0\le x \le 2)\]と幅2の長方形の畳み込みをしてみます。

import numpy as np
import matplotlib.pyplot as plt
from math import *

x = np.arange(0,2.01,0.01) 
f = x**2  //f(x)の配列を作る
g = np.ones(200)/(len(np.ones(200))) //g(x)の配列を作る
con = np.convolve(f,g) //畳み込み
conx = np.arange(-1,3,0.01)
plt.plot(conx,con)
plt.show()

まずf(x)の配列を作ります。

x = np.arange(0,2.01,0.01) 
f = x**2  //f(x)の配列を作る

ここでは0~2まで0.01刻みの配列にしています。
次にg(x)の配列を作ります。

g = np.ones(200)/(len(np.ones(200))) //g(x)の配列を作る

np.convolve関数は配列同士を片方逆順にして、かけ算して足し合わせていくというもの(([1,2,3],[0,1,0.5])の配列をnp.convolveすると([1*0,2*0+1*1,3*0+2*1+1*0.5,1*3+2*0.5,3*0.5]) = ([0,1,2.5,4,1.5])という配列ができる)なので、畳み込みにするにはどちらかを規格化しておく必要があります。
ここではg(x)を規格化して配列の和が1になるようにしています。

con = np.convolve(f,g) //畳み込み
conx = np.arange(-1,3,0.01)
plt.plot(conx,con)
plt.show()

こうして作った配列をnp.convolveで畳み込みます。
np.convolve(f,g)の配列の長さはlen(f)+len(g)-1となります。また0≦x≦2のf(x)に幅2の長方形g(x)を畳み込むので出てくる関数は-1≦x≦3です。したがってそのような配列np.arange(-1,3,0.01)を作っています。これですべて材料が揃ったので後はプロットするだけです。

出てくるグラフは次のようになります。

初めてのmatplotlib(2次関数のプロット)

pythonのグラフ作成にはmatplotlibというパッケージを利用することができます。
pythonのコードの中でグラフを作れるという点が非常に便利で、データ処理とグラフ作成を同じコード内で行うことができます。

グラフ描写にはgnuplotなんかも使えるのですがgnuplotを使う場合データ処理は別の言語で書いて、それをテキストファイルに出力して、gnuplotでそれを読み込んで描写させるという手順が必要です。matplotlibでは余計なファイルが無駄に増えないといった利点があります。またpythonのコードを使ってグラフが作れるのでpythonが使える人なら扱いやすいと思います。

早速matplotlibを使って簡単なグラフを描写してみましょう。
まずmatplotlib.pyplotとnumpyをインポートします。matplotlib.pyplotはグラフをプロットするのに必要なライブラリでpltと省略されることが多いです。またnumpyは数値計算ライブラリで行列計算したり、数値のリストを生成することができます。npと省略されることが多いです。どんなグラフを描くにしてもこの２つのライブラリはほぼ必須です。(matplotlib,numpyはインストールが必要です。インストールの仕方は他の記事を参考にしてください。)

さてmatplotlibを使って2次関数のグラフを描いてみましょう。

# -*- coding: utf-8 -*-
#y = x^2の2次関数のグラフを作成

import matplotlib.pyplot as plt  
import numpy as np 
x = np.arange(0,10,0.01) 
y = x**2
plt.plot(x,y) 
plt.show() 

図のようなグラフが表示されるはずです。
コードの説明を少しすると、最初の2行でmatplotlib.pyplotとnumpyをインストールしています。
pythonではplt.plot()関数でグラフを描写できるのですが、そこに渡すのはリストである必要があります。
今回はx = np.arange(0,10,0.01)で0から10まで0.01刻みのリストを作っています。
x = [0,0.01,0.02,...,9.99]です。その後、このリストを2乗したリストを作ってyに入れています。y = [0,0.01,0.04,...,99.8001]です。
そしてこの2つのリストをplt.plot()に渡してグラフを作っています。

リストを作るならnumpyを使わなくてもpythonのrange()関数を使えばいいと思われるかもしれません。しかし、range()関数は1以下の刻みを作ることができないので曲線の描写には不向きです。
きちんとnumpyを使ってリストを作るようにしましょう。