Posted at

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

More than 3 years have passed since last update.

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

# scipy.stats.dirichlet

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

```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()

```

# scipy.stats.wishert, scipy.stats.invwishert

こちらもMCMC等で使われることのあるウィシャート分布と逆ウィシャート分布です。

```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

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

# scipy.optimize.rosen

```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()
```

# scipy.signal.lti

```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()
```

# scipy.spatial.KDTree

```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]
```

```[(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平均法です。

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

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

# scipy.constants.g

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

# scipy.ndimage.gaussian_filter

```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()
```

# まとめ

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