Numpyの基礎（７）便利なライブラリ群
import numpy as np
import matplotlib.pyplot as plt
# jupyter notebookに画像を表示
%matplotlib inline
C言語ライブラリのラッパー¶
Numpyでは、C言語のライブラリをPythonから簡単に呼び出すことができるモジュールが用意されています。
このように、他言語のライブラリの関数をPythonで呼び出せるようにしたものをラッパー（wrapper）といいます。
もともとのライブラリをPython関数で包み込むイメージです。
モジュール名をクリックするとNumpy公式ドキュメントに飛ぶようになっています。
線形代数（np.linalg）¶
np.linalgは、元々Fortranで書かれたLAPACKがC言語に変換されたものをラップしたものです。
LAPACKの関数と一対一に対応しているわけではなく、より単純なインターフェイスで呼び出せるようになっています。
逆行列を求めたり、固有値分解、特異値分解など線形代数を網羅するライブラリです（詳細は公式ドキュメントを参照）。
# ３×３の二次元行列
A = np.random.rand(3, 3)
A
# 逆行列を求めるinv関数
B = np.linalg.inv(A)
B
# 積を取ると単位行列になる
A @ B
ある一つの軸から一つの値を取り出す操作としてノルムや最大値などがあります。
# normを取る（ようするにL2ユークリッドノルム）
np.linalg.norm([3., 4.])
このnorm関数に多次元配列を渡し、軸を指定すると、その軸に沿ってノルムを取ります。
つまりその軸がスカラーに縮約され、軸がなくなります。
# 多次元配列を定義
arr = np.arange(15).reshape(5, 3)
arr
# 各列についてノルム関数を適用する
np.linalg.norm(arr, axis=0)
# 各行についてノルム関数を適用する
np.linalg.norm(arr, axis=1)
np.linalgを用いた最小二乗法¶
下のような線形に増加するデータを最小二乗法でフィッティングすることを考えます。
t = np.linspace(0, 1, 100)
signal = 0.5 * t - 0.3 + 0.1 * np.random.randn(t.size)
plt.plot(t, signal)
plt.xlabel('Time', fontsize=20)
plt.show()
データ$d$、Green関数$G$、パラメータ$m$を使うと
$$d = Gm$$
のようになります。
これを$m$について解くと
$$m = (G^T G)^{-1}G^Td$$
# Green関数
G = np.vstack((np.ones(t.size), t)).T
G[:10]
m = np.linalg.inv(G.T @ G) @ G.T @ signal
m
plt.plot(t, signal)
plt.plot(t, m[0] + m[1] * t)
plt.xlabel('t [s]', fontsize=20)
plt.show()
うまく線形フィッティングできました。
乱数生成（np.random）¶
科学技術計算で広く使われている擬似乱数生成器であるメルセンヌ・ツイスターをラップしているモジュールです。
高品質でかつ、高速な乱数生成で有名なライブラリです。
np.randomでは、様々な確率分布で乱数を生成できるような関数群が用意されています（Distributions）。
主な関数を表にまとめました。
|関数名|説明| |:-----:|:--:| |rand|一様乱数| |randn|標準正規分布乱数| |randint|整数乱数| |choice|与えられたndarrayからランダムに選び出す| |normal|正規分布乱数|
どの関数もshapeまたはsizeが指定でき、簡単に乱数を生成できます。
# [0, 1)の範囲で一様乱数を求める
np.random.rand(10)
# ３×３の標準正規分布乱数
np.random.randn(3, 3)
# 任意のオブジェクトのリストからランダムに取り出す
np.random.choice(list("abcde"), size=10)
# ポアソン分布のヒストグラム
plt.hist(np.random.poisson(lam=[1, 4, 10], size=(10000, 3)), density=True, bins=25)
plt.show()
フーリエ変換（np.fft）¶
Cooley-Tukey型の高速離散フーリエ変換（FFT)のラッパーです。
フーリエ変換一般に言えることですが、変換と逆変換の係数が実装によって異なります（定数倍するだけですが）。
実装の詳細を知りたい人は Implementation details を参照してください。
下の表は様々な関数をまとめたものです。
|関数名|説明| |:-:|:-:| |fft|１次元FFT| |ifft|↑の逆変換| |rfft|１次元実FFT| |irfft|↑の逆変換| |fft2|２次元FFT| |ifft2|↑の逆変換| |rfft2|２次元実FFT| |irfft2|↑の逆変換|
# FFT
np.fft.fft([1., 3, 5, 3, 1])
# 実FFT
res = np.fft.rfft([1., 3, 5, 3, 1])
res
# 逆実FFT
np.fft.irfft(res)
# 逆実フーリエ変換
np.fft.irfft(res, n=5)
信号処理の一例です。
フィルターなどの高度な信号処理はscipy.signalを用いてください。
# サンプリング周波数（ナイキスト周波数＝10Hz）
rate = 20 # [1/s]
# サンプル時刻 [s]
t = np.arange(0, 10, 1/rate)
# 2Hzと6Hzにピークを持つシグナル（に正規分布ノイズを加えたもの）
signal = np.sin(2*np.pi*2*t) + np.sin(2*np.pi*6*t) + 0.1 * np.random.randn(t.size)
plt.plot(t, signal)
plt.xlim(0, 5)
plt.xlabel('time [s]', fontsize=20)
plt.title('Time series of signal', fontsize=20)
plt.show()
# power spectrum
p = np.abs(np.fft.rfft(signal))
# サンプリングレート（の逆数）を与えると各成分の周波数を返す
f = np.fft.rfftfreq(signal.size, d=1./rate)
plt.xlabel("Frequency [Hz]", fontsize=20)
plt.ylabel("Power spectrum", fontsize=20)
plt.yscale("log")
plt.plot(f, p)
plt.show()
