Scipyで見つけた数学関数あれこれ

2015年を振り返ってscipyで見つけた数学関数あれこれをまとめたいと思います。

# scipy.stats.dirichlet
LDAなどで使われるディリクレ分布の関数です。

```py
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import matplotlib.tri as tri

alphas = np.array([3.0, 4.0, 5.0])
dc = ss.dirichlet(alphas)

corners = np.array([[0.0, 0.0], [1.0, 0.0], [0.5, np.sqrt(3.0) / 2.0]])
midpoints = [(corners[(i + 1) % 3] + corners[(i + 2) % 3]) / 2.0 for i in range(3)]

def xy2bc(xy):
    s = [(corners[i] - midpoints[i]).dot(xy - midpoints[i]) / 0.75 for i in range(3)]
    return np.clip(s, 0.0, 1.0)
refiner = tri.UniformTriRefiner(tri.Triangulation(corners[:, 0], corners[:, 1]))
trimesh = refiner.refine_triangulation(subdiv=8)
pvals = [dc.pdf(xy2bc(xy)) for xy in zip(trimesh.x, trimesh.y)]
plt.tricontourf(trimesh, pvals, 200)
plt.axis('equal')
plt.show()

```

![dilichlet.png](https://qiita-image-store.s3.amazonaws.com/0/14516/81699e2c-a988-a5ce-ac61-474a8c98ddae.png)

# scipy.stats.wishert, scipy.stats.invwishert
こちらもMCMC等で使われることのあるウィシャート分布と逆ウィシャート分布です。
以下の例では、2次元の分散行列の乱数をそれぞれ10個生成します。

```py
import numpy as np
import scipy.stats as ss

w = ss.wishart(df=3, scale=np.matrix([[1.0, 0.5], [0.5, 1.0]]))
print w.rvs(10)
iw = ss.invwishart(df=3, scale=np.matrix([[1.0, 0.5], [0.5, 1.0]]))
print iw.rvs(10)
```

# scipy.special.comb, scipy.special.perm
順列と組み合わせの計算です。

```pycon
>>> import scipy.special as ss
>>> ss.comb(6, 3) # 6C3
20.0
>>> ss.perm(6, 3) # 6P3
120.0
```

# scipy.optimize.rosen
最適化のお題で使われるローゼンブロック関数です。

```py
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt

from scipy.optimize import rosen
import numpy as np

fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(np.arange(-3., 3., 0.1), np.arange(-3., 3., 0.1))
Z = [[rosen((x, y)) for x, y in zip(xx, yy)] for xx, yy in zip(X, Y)]
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
plt.show()
```
![rosen.png](https://qiita-image-store.s3.amazonaws.com/0/14516/84054401-4588-83df-3145-6a0b66420d7b.png)

# scipy.signal.lti
制御工学でお世話になる線形時不変システム(Linear Time Invariant)です。
以下の例では線形時不変システムを作成し、そのボード線図を描画しています。

```py
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# 1 / s^2 + 0.1s + 1
s = signal.lti([1], [1.0, 0.1, 1.0])
w, mag, phase = signal.bode(s, np.arange(0.1, 5.0, 0.01).tolist())

fig, (ax0, ax1) = plt.subplots(nrows=2, sharex=True)
ax0.semilogx(w, mag, 'b-')
ax0.set_ylabel("Magnitude")
ax1.semilogx(w, phase, 'r-')
ax1.set_ylabel("Phase")
plt.show()
```

![lti_bode.png](https://qiita-image-store.s3.amazonaws.com/0/14516/359d0be8-4bd2-a951-3b29-bdf691399678.png)

# scipy.spatial.KDTree
空間分割アルゴリズムのKDツリーです。
以下の例では3次元空間にグリッド状に点をばらまき、(20,20,20)の点から1.0の距離範囲にいる点をすべて抽出しています。

```py
import numpy as np
import scipy.spatial as ss                                      

x, y, z= np.mgrid[0:100, 0:100, 0:100]
points = zip(x.ravel(), y.ravel(), z.ravel())
tree = ss.KDTree(points)
a = tree.query_ball_point([20, 20, 20], 1.0)
print [points[i] for i in a]
```
```pycon
[(19, 20, 20), (20, 19, 20), (20, 20, 19), (20, 20, 20), (20, 20, 21), (20, 21, 20), (21, 20, 20)]
```

# scipy.cluster.vq.kmeans2
クラスタリングアルゴリズムのk平均法です。

```py
import numpy
from scipy.cluster.vq import kmeans2, whiten

features = numpy.array([[ 1.9,2.3],
                        [ 1.5,2.5],
                        [ 0.8,0.6],
                        [ 0.4,1.8],
                        [ 0.1,0.1],
                        [ 0.2,1.8],
                        [ 2.0,0.5],
                        [ 0.3,1.5],
                        [ 1.0,1.0]])
wf = whiten(features)  # 正規化
print kmeans2(wf, k=2)
```
中心地と各値のラベルが得られます。

```pycon
(array([[ 1.40584568,  0.69587293],
       [ 1.24002799,  2.50514254]]), array([1, 1, 0, 1, 0, 1, 0, 1, 0]))
```

# scipy.constants.g
関数ではないですが、重力加速度です。他にもいろんな物理パラメタが定義されています。

```pycon
>>> import scipy.constants
>>> scipy.constants.g
9.80665
```

# scipy.ndimage.gaussian_filter
画像処理で使われるフィルタアルゴリズムです。

```py
import scipy
from scipy import ndimage
import matplotlib.pyplot as plt

img = scipy.misc.lena().astype(float)
fimg = ndimage.gaussian_filter(img, 3)

plt.figure(figsize=(8, 4))
plt.subplot(121)
plt.imshow(img, cmap=plt.cm.gray)
plt.subplot(122)
plt.imshow(fimg, cmap=plt.cm.gray)
plt.show()
```

![lena.png](https://qiita-image-store.s3.amazonaws.com/0/14516/d50626c2-ecfb-4700-66c6-a4fb25d642d0.png)

# まとめ
scipyはいろんな分野のいろんな関数やアルゴリズムが実装されていて、リファレンスを眺めているだけでも楽しいですね。来年もお世話になります。