科学技術計算のためのパッケージ

SciPy http://www.scipy.org/ は様々な科学技術計算のための実装を含むパッケージ(モジュール群)である.

  • optimize: 非線形最適化と方程式の根
  • spatial: 計算幾何学
  • stats: 確率・統計
  • interpolate: 補間
  • integrate: 積分

最適化


SciPyの最適化サブモジュールoptimize

非線形最適化ならびに非線形方程式の根を求めるための解法群

optimizeに含まれる非線形最適化のための関数

1 変数用のminimize_scalar

例: 1変数関数の最小化

関数 $f(x)= x/(1+x^2)$

$x=-1$ のとき最小値 $-1/2$ をとる.

  • グラフによる確認

  • Brent法

  • 黄金分割法

  • 引数boundsで範囲を与えたBrent法

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt


def f(x):
    return x / (1 + x**2)


x = np.linspace(-10, 10, 100)
y = f(x)
plt.plot(x, y)
from scipy.optimize import minimize_scalar

print("Brent=========================")
print(minimize_scalar(f, method="Brent"))

print("Golden=========================")
print(minimize_scalar(f, method="Golden"))

print("Bounded=========================")
print(minimize_scalar(f, bounds=(-2.0, -1.5), method="Bounded"))
Brent=========================
     fun: -0.5
    nfev: 14
     nit: 9
 success: True
       x: -0.9999999915247861
Golden=========================
     fun: -0.5
    nfev: 45
     nit: 39
 success: True
       x: -0.9999999904634918
Bounded=========================
     fun: -0.4615379352317438
 message: 'Solution found.'
    nfev: 24
  status: 0
 success: True
       x: -1.500004447286288

minimize(f,x0)


複数の変数をもつ関数fを初期解x0からの探索で最小化

引数methodで探索のためのアルゴリズムを設定

SciPyのoptimizeサブモジュールの minimize(f,x0)


複数の変数をもつ関数fを初期解x0からの探索で最小化

引数methodで探索のためのアルゴリズムを設定

  • 'Nelder-Mead':Nelder–Mead法(単体法)

  • 'Powell':Powellの共役方向法

  • 'CG':共役勾配法(conjugate gradient)

  • 'BFGS':Broyden–Fletcher–Goldfarb–Shanno(BFGS)法

  • 'Newton-CG':Newton共役勾配法

  • 'L-BFGS-B':記憶制限付きBFGS法

  • 'TNC':打ち切りNewton共約勾配法

  • 'COBYLA':Constrained Optimization BY Linear Approximation

  • 'SLSQP':Sequential Least SQuares Programming

  • 'dogleg':ドッグレッグ信頼領域法

  • 'trust-ncg':信頼領域Newton共約勾配法

他の引数

  • jac : 勾配ベクトル(Jacobian)

  • hess : Hesse行列(Hessian)

  • bounds :上下限(限界値)制約

  • constraints: 一般の制約式($=0$ もしくは $\geq 0$)

    • 制約の種類は文字列 'type' をキーとし,値を制約の種類を表す文字列 ('eq', 'ineq').

    • 制約関数は文字列 'fun'をキーとし,制約を表す関数を値

    • Jacobianは文字列 'jac'をキーとし,Jacobianを表す関数を値

制約の例:不等式 $x_0 +2x_1 \geq 10$ を表す制約

{'type': 'ineq', 'fun': lambda x: x[0] + 2*x[1] -10}

制約なし最適化の例

Rosenbrock関数

$$f(x_1, x_2, \ldots, x_n) = \sum_{i=1}^{n-1} (100(x_i^2 - x_{i+1})^2 + (1-x_i)^2)$$


  • 関数: rosen
  • Jacovian: rosen_der
  • Hessian: rosen_hess
    sum(100.0*(x[1:] - x[:-1]**2.0)**2.0 + (1 - x[:-1])**2.0)

x=(1,1,...,1) のとき最小

  • 2次元の場合の描画
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import LogNorm
import matplotlib.cm as cm

ax = Axes3D(plt.figure())
s = 0.05
X = np.arange(-2, 2.0 + s, s)
Y = np.arange(-1, 3.0 + s, s)
X, Y = np.meshgrid(X, Y)
Z = (1.0 - X) ** 2 + 100.0 * (Y - X * X) ** 2
ax.plot_surface(
    X, Y, Z, rstride=1, cstride=1, norm=LogNorm(), cmap=cm.jet, linewidth=0.2
)
None
/var/folders/69/5y96sdc94jxf6khgc8mlmxrr0000gn/T/ipykernel_73092/2668126910.py:5: MatplotlibDeprecationWarning: Axes3D(fig) adding itself to the figure is deprecated since 3.4. Pass the keyword argument auto_add_to_figure=False and use fig.add_axes(ax) to suppress this warning. The default value of auto_add_to_figure will change to False in mpl3.5 and True values will no longer work in 3.6.  This is consistent with other Axes classes.
  ax = Axes3D(plt.figure())
import plotly.graph_objs as go
import plotly

X, Y = np.mgrid[-2:2:0.05, -1:3:0.05]
Z = (1.0 - X) ** 2 + 100.0 * (Y - X * X) ** 2
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=np.log(Z + 0.01)))
# fig.show()
# plotly.offline.plot(fig); #Jupyter Labの場合

5次元の場合の比較

  • 単体法(Nelder-Mead法:勾配情報なし)

  • BFGS(準Newton法)(勾配情報のみ使う)

  • 信頼領域Newton共役勾配法(勾配とHessian)

from scipy.optimize import minimize, rosen, rosen_der, rosen_hess

x0 = [2.0 for i in range(5)]
res = minimize(rosen, x0, method="Nelder-Mead")
print(res.x)

res = minimize(rosen, x0, jac=rosen_der, method="BFGS")
print(res.x)

res = minimize(rosen, x0, jac=rosen_der, hess=rosen_hess, method="trust-ncg")
print(res.x)
[1.00000047 1.00000079 1.0000027  1.00000292 1.00000473]
[0.99999999 0.99999998 0.99999995 0.99999992 0.99999985]
[1.         1.         1.         0.99999999 0.99999993]

Beale 関数


$$ f(x_1, x_2) = \displaystyle\sum_{i=1}^{3} (c_i - x_1 (1- (x_2)^i )^2 $$

$c_1=1.5, c_2=2.25, c_3=2.625$

$x_2=1$ 付近で湾曲

$x=(3,0.5)$ のとき最小値 $f^* =0$

  • 単体法(Nelder-Mead法)

  • BFGS(準Newton法)

ax = Axes3D(plt.figure())
c = {1: 1.5, 2: 2.25, 3: 2.625}
s = 0.5
X = np.arange(-5, 5.0 + s, s)
Y = np.arange(-5, 5.0 + s, s)
X, Y = np.meshgrid(X, Y)
Z = sum((c[i] - X * (1 - Y**i)) ** 2 for i in range(1, 4))
ax.plot_surface(
    X, Y, Z, rstride=1, cstride=1, norm=LogNorm(), cmap=cm.jet, linewidth=0.2
)
/var/folders/69/5y96sdc94jxf6khgc8mlmxrr0000gn/T/ipykernel_73092/1353534371.py:1: MatplotlibDeprecationWarning:

Axes3D(fig) adding itself to the figure is deprecated since 3.4. Pass the keyword argument auto_add_to_figure=False and use fig.add_axes(ax) to suppress this warning. The default value of auto_add_to_figure will change to False in mpl3.5 and True values will no longer work in 3.6.  This is consistent with other Axes classes.

c = {1: 1.5, 2: 2.25, 3: 2.625}
X, Y = np.mgrid[-5:5:0.05, -5:5:0.05]
Z = sum((c[i] - X * (1 - Y**i)) ** 2 for i in range(1, 4))
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=np.log(Z + 0.01)))
# fig.show();
# plotly.offline.plot(fig);  #Jupyter Labの場合

Powell関数


$$ f(x_1, x_2, x_3, x_4) = (x_1-10x_2)^2 + 5(x_3-x_4)^2 + (x_2-2x_3)^4+ 10(x_1-x_4)^4$$

最適解 $(0,0,0,0)$, $f^*=0$

最適解におけるHesse行列が特異

Ackley関数


$$f(x,y) = -20\exp\left(-0.2\sqrt{0.5\left(x^{2}+y^{2}\right)}\right)-\exp\left(0.5\left(\cos\left(2\pi x\right)+\cos\left(2\pi y\right)\right)\right) + e + 20$$

最適解 $x=(0,0), f^*= 0$

$-5 \le x,y \le 5$

ax = Axes3D(plt.figure())
s = 0.1
X = np.arange(-3, 3.0 + s, s)
Y = np.arange(-3, 3.0 + s, s)
X, Y = np.meshgrid(X, Y)
Z = (
    -20 * np.exp(-0.2 * np.sqrt(0.5 * (X**2 + Y**2)))
    - np.exp(0.5 * (np.cos(2 * np.pi * X) + np.cos(2 * np.pi * Y)))
    + np.exp(1)
    + 20
)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet, linewidth=0.2)
None
/var/folders/69/5y96sdc94jxf6khgc8mlmxrr0000gn/T/ipykernel_73092/3623169322.py:1: MatplotlibDeprecationWarning:

Axes3D(fig) adding itself to the figure is deprecated since 3.4. Pass the keyword argument auto_add_to_figure=False and use fig.add_axes(ax) to suppress this warning. The default value of auto_add_to_figure will change to False in mpl3.5 and True values will no longer work in 3.6.  This is consistent with other Axes classes.

X, Y = np.mgrid[-3:3:0.05, -3:3:0.05]
Z = (
    -20 * np.exp(-0.2 * np.sqrt(0.5 * (X**2 + Y**2)))
    - np.exp(0.5 * (np.cos(2 * np.pi * X) + np.cos(2 * np.pi * Y)))
    + np.exp(1)
    + 20
)
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=Z))
# fig.show()
# plotly.offline.plot(fig); #Jupyter Labの場合

問題

$f(x,y) = x^2 +xy+y^2$ を最小にする解を共役勾配法(引数methodは'CG')を用いて最小化せよ.

X, Y = np.mgrid[-3:3:0.05, -3:3:0.05]
Z = X**2 + X * Y + Y**2
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=Z))
# fig.show()
# plotly.offline.plot(fig);  #Jupyter Labの場合

問題

$f(x,y)=\sqrt{1+ x^2+y^2}$をNelder–Mead法(引数methodは'Nelder-Mead')とBFGS法(引数methodは'BFGS')でを用いて最小化せよ. 解の精度を比較せよ.

X, Y = np.mgrid[-3:3:0.05, -3:3:0.05]
Z = np.sqrt(1 + X**2 + Y**2)
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=Z))
# fig.show()
# plotly.offline.plot(fig);  #Jupyter Labの場合

問題

$f(x,y) = x^3 -3xy + y^3$ をBFGS法(引数methodは'BFGS')で最小化せよ. 初期点を色々変えてみて,どのような解が得られるか観察せよ.

X, Y = np.mgrid[-3:3:0.05, -3:3:0.05]
Z = X**3 - 3 * X * Y + Y**3
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=Z))
# fig.show()
# plotly.offline.plot(fig);  #Jupyter Labの場合

問題(動的価格付け問題)

あなたはスーパーの店長だ.普段100円で販売している卵の値段を変えて収益を最大化することを考えている. 現在の価格での販売量は1日50個で,1円値下げすると2個余計に売れるようになり,逆に1円値上げすると2個売れないようになると推測されている. さて,おおよそ何円で売れば良いだろうか?

問題(経済発注量問題)

あなたは工場の在庫管理部長だ.いま,ある品目の発注費用が1回あたり1000円で,在庫費用が1日あたり1個あたり1円と推定されている. 1日あたりの需要量が100個としたとき,何日おきに発注すれば良いだろうか?

変数は品目の発注間隔(サイクル時間)であり,それを $x$ とすると,1日あたりの発注費用は $\frac{1000}{x}$, 在庫費用は $\frac{1 \cdot 100 \cdot x}{2}$ と書くことができる.

発注費用と在庫費用の和を最小化する発注間隔を求めよ.ただし,発注間隔は整数でなくても良いものとする.

問題(放物線)

鳥が放物線 $y=x^2 + 10$ を描いて飛んでいる.いま $(10,0)$ の位置にいるカメラマンが,鳥との距離が最小の地点で シャッターを押そうとしている.鳥がどの座標に来たときにシャッターを押せば良いだろうか?

問題(楕円軌道)

彗星が地球に接近している.地球を原点 $(0,0,0)$ としたとき彗星の描く軌道は曲面 $2x^2 + y^2 + z^2 = 1000$ 上にあることが 予測されている.地球に最も接近するときの距離を求めよ.

問題(体積)

断面の面積の和が $100$ cm${}^2$ の直方体で,体積最大のものは何か?また体積最小のものは何か?

制約付き最適化の例


半径 $\sqrt{2}$ の円内の制約の下で $x+y$ を最大化

  • COBYLA (線形近似法;Constrained Optimization BY Linear Approximation)で求解

  • SLSQP(逐次最小2乗法;Sequential Least SQuares Programming)で求解

目的関数や制約式は,普通の関数(defで定義)で書いても,lambda関数として書いても良い.以下の例では2通りの書き方で目的関数を書いている.(SciPyでの目的関数は最小化なので,-1を乗じてある.)

制約は辞書を用いて定義する. 辞書には,以下のキーと値の対を入れる.

  • 'type': 制約のタイプで等式のときは 'eq',不等式のときには 'ineq'を入れる.
  • 'fun': 制約を表す関数を入れる. 制約は必ず,$0$ 以上の制約に直してから入れる.

以下の例では, $x^2 + y^2 \leq 2$ を $-x^2 -y^2 +2 \geq 0$ に直してから入れている.

その後に,minimmize関数で最小化するが,制約が複数ある場合には,制約(辞書)のリストを constraints引数に入れる.

from scipy.optimize import minimize


def f(x):
    return -x[0] - x[1]


bnds = ((0, None), (0, None))  # 上下限制約の入れ方
f = lambda x: -x[0] - x[1]
con = {"type": "ineq", "fun": lambda x: -x[0] ** 2 - x[1] ** 2 + 2}
res = minimize(fun=f, x0=[2, 2], method="COBYLA", constraints=con)
print("COBYLA========= \n", res)
res = minimize(fun=f, x0=[2, 2], method="SLSQP", constraints=con)
print("SLSQP========== \n", res)
COBYLA========= 
      fun: -2.000000006134913
   maxcv: 2.0000961775679116e-08
 message: 'Optimization terminated successfully.'
    nfev: 47
  status: 1
 success: True
       x: array([1.00006218, 0.99993783])
SLSQP========== 
      fun: -2.000000092931198
     jac: array([-1., -1.])
 message: 'Optimization terminated successfully'
    nfev: 15
     nit: 5
    njev: 5
  status: 0
 success: True
       x: array([1.00000005, 1.00000005])

問題(制約付きWeber問題)

以下のような位置と人数をもった7件の家がある.

x座標 y座標 人数
1 24 54 2
2 60 63 1
3 1 84 2
4 23 100 3
5 84 48 4
6 15 64 5
7 52 74 4

みんなで出資して新しく井戸を掘ることになったが, 話し合いの末 「各家が水を運ぶ距離 $\times$ 各家の水の消費量」の総和が最小になる場所に 井戸を掘ることにした. ただし, 各家の水の消費量は人数に比例するものとする.

  1. どこを掘っても水が出るものとしたとき,どのようにして掘る場所を決めれば良いだろうか?

  2. $(50,50)$ から半径 $40$ の円状の領域が砂漠になっていて,そこには井戸が掘れないときにはどうしたら良いだろうか?

家の集合を $H$,家 $i (\in H)$ の位置を$(x_i, y_i)$ とする. 井戸の位置を $(X, Y)$ とすれば,家 $i$ から井戸までの距離は $$ \sqrt{(x_i - X)^2 + (y_i - Y)^2} $$ である.家 $i$ が $1$日に必要とする水の量(人数)を $w_i$ としたとき, $$ \sum_{i \in H} w_i \sqrt{(x_i-X)^2 + (y_i-Y)^2} $$ を最小にする $(X, Y)$ を求めよ.

定式化: $$ \begin{array}{ll} minimize & \sum_{i \in H} w_i \sqrt{(x_i-X)^2 + (y_i-Y)^2} \\ s.t. & (50-X)^2 +(50-Y)^2 \geq 40^2 \end{array} $$

x = [24, 60, 1, 23, 84, 15, 62]
y = [54, 63, 84, 100, 48, 64, 74]
w = [2, 1, 2, 3, 4, 5, 4]

問題(ポートフォリオ最適化)

100万円のお金を5つの株に分散投資したいと考えている. 株の価格は,現在はすべて $1$株あたり $1$円だが, 証券アナリストの報告によると,それらの株の $1$年後の価格と標準偏差はそれぞれ 以下の表のように確率的に変動すると予測されている.

株式 1 2 3 4 5
期待値 ($r_i$) 1.01 1.05 1.08 1.10 1.20
標準偏差 ($\sigma$) 0.07 0.09 0.1 0.2 0.3

目的は$1$年後の資産価値を最大化することである. しかしながら,よく知られているように,1つの株式に集中投資するのは 危険であり,大損をすることがある. 期待利回りを $\alpha$%以上としたとき,標準偏差の2乗和を最小にするように投資するにはどうすれば良いだろうか?

ポートフォリオ最適化(定式化):

$$ \begin{array}{lrcl} minimize & \sum_{i=1}^n \sigma_i^2 x_i^2 \\ s.t. & \sum_{i=1}^n r_i x_i &\geq& \alpha \\ & \sum_{i=1}^n x_i &=& 1 \\ & x_i &\geq& 0 \ \ \ \forall i=1, 2, \ldots, n \end{array} $$

問題(交通量割り当て問題)

  • 4 本の道への交通量配分
  • 5000台の車が移動する総時間を最小化

$移動時間 = 基本移動時間 \times \left(1+ \left(\frac{台数}{交通容量}\right)^4 \right)$

基本移動時間 交通容量
1 15 1000
2 20 2000
3 30 3000
4 35 4000

問題(学区割り当て問題) (難)

  • 2 つの学校(定員100人)を作る
  • 学生たちの歩く距離の合計を最小
  • 新しい学校の位置と割当
学区 x座標 $x_i$ y座標 $y_i$ 人数 $n_i$
1 0 0 40
2 0 100 40
3 100 0 40
4 100 100 40
5 50 50 40

学区割り当て問題(定式化)

$$ \begin{array}{l l l} minimize & \sum_{i \in I} \sum_{j=0}^1 z_{ij} \sqrt{ (x_i-X_j)^2+(y_i-Y_j)^2 } & \\ s.t. & \sum_{i \in I} z_{ij} \leq 100 & \forall j=0,1 \\ & \sum_{j=0}^1 z_{ij} = n_i & \forall i \in I \\ & z_{ij} \geq 0 & \forall i \in I; j=0,1 \end{array} $$
  • 非凸関数の最適化問題
  • 適切な初期解を与えないと妥当な解が得られない可能性

問題(複数品目の経済発注量問題)

古典的な経済発注量モデルの複数品目への拡張(容量制約付き)

  • $I$ : 品目の集合
  • $w_i$ : 品目 $i~(\in I)$ の大きさ
  • $W$ : 倉庫に置ける量の上限.
  • $F_i$ : 品目 $i$ の発注費用
  • $h_i$ : 品目 $i$ の在庫保管費用
  • $d_i$ : 品目 $i$ 需要量

容量制約を破らない条件の下で,総費用が最小になる発注方策を求める.

複数品目の経済発注量問題(定式化)

$T_i$:品目 $i$ の発注間隔(発注頻度の逆数)

$$ \begin{array}{l l l} minimize &\sum_{i \in I} \frac{F_i}{T_i} + \frac{h_i d_i T_i}{2} & \\ s.t. & \sum_{i \in I} w_i d_i T_i \leq W & \\ & T_i > 0 & \forall i \in I \end{array} $$

数値は自分で適当に定めて解け.

曲線へのあてはめ(最小2乗法)

  • curve_fit(f,x,y): 配列 $x,y$ で与えたデータに対して誤差の最小2 乗和を最小にするような関数 $f$ のパラメータを求める

    返値:

    • 最小2 乗和を最小にするような推定パラメータを表す配列

    • 共分散行列を表す2 次元配列のタプル

例: 2 次関数 $a x^2+ bx+c$ へのあてはめ

from scipy.optimize import curve_fit


def f(x, a, b, c):
    return a * x**2 + b * x + c


xdata = np.linspace(-10, 10, 20)
ydata = f(xdata, 1, 1, 1) + np.random.normal(size=len(xdata))
param, cov = curve_fit(f, xdata, ydata)
print(param)
[1.00359825 0.98692547 1.03629323]
plt.plot(xdata, ydata, "ro")
plt.plot(xdata, f(xdata, param[0], param[1], param[2]))
None

問題(コーヒーショップの売上)

あなたはコーヒーショップを開店した.オープン当初の5日間の売り上げは, $15,20,34,42,58$ と右肩上がりである.1日あたりどのくらい増えているかを 線形関数 $y=ax+b$ へ当てはめることによって予測せよ.

1変数関数の根

$x^2 -4 \sin x = 0$ の解($x^2 -4 \sin x$ の根)を求める.

$x=0$ だけでなく $x=1.9$ 付近でも根をもつ.

scipy.optimizeのrootもしくはbrentq関数を用いる.

$f(x)=0$ の解を出すためには, rootは引数として関数fと初期値 x0を与え, brentqは関数fと範囲 aとbを与える.

def f(x):
    return x**2 - 4 * np.sin(x)


x = np.linspace(-1, 3, 100)
plt.plot(x, np.linspace(0, 0, 100))  # x=0
y = f(x)
plt.plot(x, y)
None
from scipy.optimize import root, brentq

print("Brentq=", brentq(f, 1.0, 3.0))
print("Root finding ======\n ", root(f, 2.5))
Brentq= 1.9337537628270214
Root finding ======
      fjac: array([[-1.]])
     fun: array([1.84741111e-13])
 message: 'The solution converged.'
    nfev: 8
     qtf: array([-3.12600879e-08])
       r: array([-5.2877011])
  status: 1
 success: True
       x: array([1.93375376])

問題(4次方程式の根)

$x^4 + x^3 -7x^2-x+6 =0$ の解をすべて求めよ.

解と図は以下のサイトで x^4 + x^3 -7x^2-x+6 =0 と入力することによって確認できる.

https://www.wolframalpha.com/

計算幾何学


計算幾何とデータ構造に関するサブモジュール spatial

計算幾何学(computational geometry)とは,幾何的問題を効率よく解くための基本アルゴリズムの体系化をめざす学問分野である. 1970年代の中頃に生まれたこの分野は, その有用性のためその後急速に発展している.

主なクラス

  • K次元木 KDTree

  • 凸包 ConvexHull

  • Delaunay 3 角形分割 Delaunay

  • Voronoi図 Voronoi

  • 距離計算のためのサブモジュール distance

$K$-d木

KDTreeは$K$次元の点データに対する$K$-d木とよばれるデータ構造のクラスである. ここで$K$-d木($K$-dimensional tree)とは, $K$次元のEuclid空間にある点を保管するための空間分割データ構造である. $K$-d木を用いることによって,範囲探索や最近点探索などを高速に行うことができる.

$K$-d木は空間の情報を2分木として保管している. ここで2分木とは,(root)とよばれる点(これは与えた点全体を表す)を2つの子点(child node) に分割する操作を繰り返すことによって生成されるデータ構造である.

2次元を例にすると,$K$-d木は以下のように構築される. 点集合を含む領域(長方形)を $x,y$軸のいずれか長い方の中央値で分割することによって2つの長方形を生成し,これらを子点とする. 子点には長方形に含まれる点を保持する. 各子点に対して同様の操作を繰り返し,長方形に含まれる点の数が一定数以下になったら,終了する. 子点をもたない点を葉(leaf)とよぶ.

コンストラクタの引数:

  • points: NumPyの形状 (npoints, ndim) をもつ浮動小数点数の配列; npointsは点数,ndimは次元数.

  • leafsize: 葉(子をもたない点)に含まれる点数を表す自然数;

メソッドquery(x)は,点(もしくは配列)xに含まれる各点からの最近点探索を行う. ここでx(の各要素)は,$K$-d木と同じ次元をもつ点の座標を表す配列とする. 返値は点xから近い順に並べ替えられた距離(の配列)と近い点の位置(を表す配列)である. オプションで与えることができる引数は以下の通り.

  • k は近い点の数を表すパラメータであり,自然数を与える.既定値は $1$ である.
  • p は距離を計算する際のノルムを表すパラメータである. $n$ 次元の点 $x,y$ 間のMinkowski $p$-ノルムは以下のように計算される. $$ \left( \sum_{i=1}^n |x_i -y_i|^p \right)^{1/p} $$

pは $1$ 以上の浮動小数点数を与え,$1$ のときは $L_1$(マンハッタン)ノルム,$2$ のときはEuclidノルム, 無限大のときは $L_{\infty}$ノルムになる.既定値は $2$ である.

例として, $[0,1]^2$ にランダムに発生させた $1000$個の2次元の点に対して, $(0.1,0.1)$ と $(0.5,0.5)$ から近い順に10個の点を表示する.

%matplotlib inline
from scipy.spatial import KDTree
import numpy as np
import matplotlib.pyplot as plt

points = np.random.rand(10000, 2)
tree = KDTree(points)
dis, near = tree.query(x=[[0.1, 0.1], [0.5, 0.5]], k=30)
print(dis)  # 各点から近い点への距離(近い順)
plt.plot(points[:, 0], points[:, 1], "b+")
plt.plot(points[near, 0], points[near, 1], "ro")
None
[[0.00233924 0.00594753 0.00723594 0.00854518 0.00940984 0.01088481
  0.01096928 0.01420615 0.01430367 0.01497642 0.01521166 0.01941017
  0.01993426 0.02001961 0.02203703 0.0222043  0.02221201 0.02414512
  0.02418772 0.02427058 0.0265167  0.02677313 0.02703475 0.0276239
  0.02823189 0.02826491 0.02858217 0.02923699 0.02940296 0.0295949 ]
 [0.00137277 0.00336361 0.00764734 0.01013536 0.01292676 0.01298923
  0.01514037 0.01586783 0.01597838 0.01693038 0.01926416 0.02008635
  0.02118614 0.02122255 0.02124094 0.02155826 0.02295749 0.02309385
  0.02336436 0.02344152 0.02396423 0.02440003 0.02500508 0.02514418
  0.02538147 0.02646367 0.02716226 0.02774051 0.02810013 0.02863846]]

メソッドquery_ball_point(x, r)は点(もしくは配列)xから距離r以内のすべての点のリスト(もしくはリストの配列)を返す.

その他の引数はqueryと同じである.

$[0,1]^2$ にランダムに発生させた $1000$個の2次元の点に対して, 中心 $(0.5,0.5)$ から $L_1$ノルムでの距離 $0.3$ 以内の点を表示を行う.

points = np.random.rand(1000, 2)
tree = KDTree(points)
ball = tree.query_ball_point(x=[0.5, 0.5], r=0.3, p=1)
plt.plot(points[:, 0], points[:, 1], "b+")
plt.plot(points[ball, 0], points[ball, 1], "ro")
None

凸包クラス ConvexHull


  • 凸包:計算幾何学の基本となる概念;与えられたEuclid空間内の点を含む最小の凸多面体

コンストラクタの引数:

  • points: NumPyの形状 (npoints, ndim) をもつ浮動小数点数の配列; npointsは点数,ndimは次元数.

凸包クラスConvexHullは以下の属性とメソッドをもつ.

  • pointsは入力された点の座標を表す配列である.
  • verticesは凸包を構成する点を表す配列である.2次元の場合には,反時計回りの順に保管されるが,3次元以上の場合には任意の順になる.
  • simplicesは凸包を構成する面(単体)を保持する配列である. 面の数の配列であり,配列の各要素は点の次元の長さの配列である.
  • neighborsは凸包を構成する面(単体)の隣接関係を表す配列である.配列の各要素は隣接する面の番号の配列である.

また,2次元の凸包を描画するための関数convex_hull_plot_2dも準備されている. これは凸包インスタンスを引数とし,matplotlibの図を返す.

from scipy.spatial import ConvexHull, convex_hull_plot_2d

points = np.random.rand(300, 2)
hull = ConvexHull(points)
convex_hull_plot_2d(hull)
print(hull.vertices)
[194  26  90 150  68  70 287 164  73 274 236 210 198 157 253]

Delaunay 3 角形分割クラス Delaunay


  • Delaunay3 角形分割: (2次元のときは)領域を点を通る3 角形に分割したとき,どの3角形の頂点を通る円も,他の点を含まない3角形分割

  • 近い点同士を見つけるのに便利;Voronoi図と双対の関係

コンストラクタの引数:

  • points: NumPyの形状 (npoints, ndim) をもつ浮動小数点数の配列; npointsは点数,ndimは次元数.

Delaunay3角形分割クラスDelaunayは以下の属性とメソッドをもつ.

  • pointsは入力された点の座標を表す配列である.
  • simplicesはDelaunay3角形分割に含まれる単体(2次元のときは3角形)を保持する配列である. 単体の数の配列であり,各要素は次元 $+1$ の長さの配列に単体の頂点番号が保管される.
  • neighborsはDelaunay3角形分割に含まれる単体の隣接関係を表す配列である.配列の各要素は, 単体内の各点の反対側に隣接する単体の番号の配列である.対応する単体がない場合には $-1$ が保管される.

また,2次元のDelaunay3角形分割を描画するための関数delaunay_plot_2dも準備されている. これはDelaunay3角形分割インスタンスを引数とし,matplotlibの図を返す.

from scipy.spatial import Delaunay, delaunay_plot_2d

np.random.seed(5)
points = np.random.rand(80, 2)
tri = Delaunay(points, furthest_site=False)
fig = delaunay_plot_2d(tri)

Voronoi図クラス Voronoi


  • Voronoi図: 各点に近い空間で領域分けされた図

  • 点に近い領域を見つけるのに便利;Delaunay 3 角形分割の双対

  • 与えられた点を母点,各母点に近い空間から成る領域をVoronoi領域. Voronoi領域の境界をVoronoi境界,Voronoi境界の交点をVoronoi点

コンストラクタの引数:

  • points: NumPyの形状 (npoints, ndim) をもつ浮動小数点数の配列; npointsは点数,ndimは次元数.

Voronoi図クラスVoronoiは以下の属性とメソッドをもつ.

  • pointsは入力された点の座標を表す配列である.
  • verticesはVoronoi点の座標を表す配列である.
  • ridge_pointsはVoronoi境界を垂直に横切る母点対を表す. 実は,これはDelaunay3角形分割の枝(辺)に対応し,これがVoronoi図とDelaunay3角形分割が互いに双対といわれる所以である.
  • ridge_verticesはVoronoi境界を表す点対である. Voronoi点の番号,もしくは点がVoronoi図の外側にある場合には $-1$ を保管する.
  • regionsはVoronoi領域を表す点の配列である.点がVoronoi図の外側にあるときには $-1$ を保管する.
  • point_regionは各母点が含まれるVoronoi領域の番号を表す.

また,2次元のVoronoi図を描画するための関数voronoi_plot_2dも準備されている. これはVoronoi図インスタンスを引数とし,matplotlibの図を返す.

例として格子上に配置した点に対するVoronoi図を求め,それを描画する.

同時に,対応するDelaunay3角形分割も示す.

from scipy.spatial import Voronoi, voronoi_plot_2d

points = np.array(
    [[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]]
)
vor = Voronoi(points)
fig = voronoi_plot_2d(vor)
print("vertices=", vor.vertices)
print("reige_points", vor.ridge_points)
vertices= [[0.5 0.5]
 [0.5 1.5]
 [1.5 0.5]
 [1.5 1.5]]
reige_points [[0 3]
 [0 1]
 [2 5]
 [2 1]
 [1 4]
 [7 8]
 [7 6]
 [7 4]
 [8 5]
 [6 3]
 [4 5]
 [4 3]]
tri = Delaunay(points)
fig = delaunay_plot_2d(tri)

Delaunay3角形分割の例題に対するVoronoi図を計算し,対応するDelaunay3角形分割と比較する.

np.random.seed(5)
points = np.random.rand(8, 2)
vor = Voronoi(points)
fig = voronoi_plot_2d(vor)

tri = Delaunay(points) fig = delaunay_plot_2d(tri)

距離計算のためのサブモジュールdistance


  • euclidean(u,v): uとvの間のEuclid距離(直線距離)

    $\|u-v\|_2 = \left(\sum_{i} (u_i-v_i)^2 \right)^{1/2}$

  • chebyshev(u,v): uとvの間のChebyshev距離

    $(L_{\infty}$ ノルム) $\max_{i} |u_i -v_i|$

  • cityblock: uとvの間のManhattan距離($L_{1}$ ノルム)

    $\sum_{i} |u_i-v_i|$

  • minkowski(u,v,p): uとvの間のMinkowski距離((Minkowski p-ノルム)

    $\| u-v \|_p =\left(\sum_{i} |u_i-v_i|^p \right)^{1/p}$

例題:最近点

  • 2 次元Euclid平面上にランダムに分布した10000個の点に対する最も近い点を求める
  • 最も近い点は,点を母点としたVoronoi図の隣接する領域の母点であることを利用
  • $K$-d木と比較
from scipy.spatial.distance import euclidean


def NN_KD(points):
    tree = KDTree(points)
    near = {}
    total_cost = 0.0
    for i, p in enumerate(points):
        dis, near_list = tree.query(p, k=2)
        near[i] = near_list[1]


def NN_Voronoi(points):
    vor = Voronoi(points)
    dis, near = {}, {}
    for i in range(len(points)):
        dis[i] = np.inf
    for (i, j) in vor.ridge_points:
        d = euclidean(points[i], points[j])
        if d < dis[i]:
            dis[i] = d
            near[i] = j
        if d < dis[j]:
            dis[j] = d
            near[j] = i


n = 10000
points = np.random.rand(n, 2)
%time NN_KD(points)
%time NN_Voronoi(points)
CPU times: user 1.15 s, sys: 3.14 ms, total: 1.16 s
Wall time: 1.16 s
CPU times: user 536 ms, sys: 3.98 ms, total: 540 ms
Wall time: 540 ms

例題:Euclid最小木問題


  • 点間の距離はEuclid距離

  • Euclid最小木問題の最適解に含まれる枝は,必ずVoronoi図の隣接する領域間の母点対の枝集合に含まれる

  • Voronoi図の隣接する領域間の母点対は ridge_points を用いて列挙

import networkx as nx


def MST(points):
    G = nx.Graph()
    for i, p in enumerate(points):
        for j, q in enumerate(points):
            if i < j:
                G.add_edge(i, j, weight=euclidean(p, q))
    E = nx.minimum_spanning_edges(G)


def MST_Voronoi(points):
    G = nx.Graph()
    vor = Voronoi(points)
    for (i, j) in vor.ridge_points:
        G.add_edge(i, j, weight=euclidean(points[i], points[j]))
    E = nx.minimum_spanning_edges(G)
    return G, E


n = 200
points = np.random.rand(n, 2)

%time MST(points)
%time MST_Voronoi(points)

G, E = MST_Voronoi(points)
pos = {i: (points[i][0], points[i][1]) for i in range(n)}
nx.draw(G, pos=pos, node_size=10, edgelist=list(E))
CPU times: user 177 ms, sys: 1.58 ms, total: 178 ms
Wall time: 178 ms
CPU times: user 7.31 ms, sys: 83 µs, total: 7.39 ms
Wall time: 7.08 ms

問題: 色々な距離

2次元座標空間上の点[0,0] と点 [1,1] の間の距離を求めよ.

  1. Euclid距離(直線距離)

  2. Chebyshev距離 $(L_{\infty}$ ノルム)

  3. Manhattan距離($L_{1}$ ノルム)

  4. $p=3$ のMinkowski距離((Minkowski p-ノルム)

問題:$K$-d木

日本の都道府県の県庁所在地の経度・緯度のリストが以下のように与えられている. これを2次元平面上の座標と考えて$K$-d木を構築し,仙台市 (140.87194,38.26889) に直線距離で近い点の番号を順に3つ求めよ.

lon = [
    141.34694,
    140.74,
    141.1525,
    140.87194,
    140.1025,
    140.36333,
    140.46778,
    140.44666999999998,
    139.88361,
    139.06083,
    139.64889,
    140.12333,
    139.69167,
    139.6425,
    139.02361000000002,
    137.21139,
    136.62556,
    136.22194,
    138.56833,
    138.18111000000002,
    136.72222,
    138.38306,
    136.90667,
    136.50861,
    135.86833,
    135.75556,
    135.52,
    135.18306,
    135.83278,
    135.1675,
    134.23833,
    133.05056000000002,
    133.935,
    132.45944,
    131.47139,
    134.55944,
    134.04333,
    132.76611,
    133.53111,
    130.41806,
    130.29889,
    129.87361,
    130.74167,
    131.6125,
    131.42389,
    130.55806,
    127.68111,
]
lat = [
    43.06417,
    40.82444,
    39.70361,
    38.26889,
    39.71861,
    38.24056,
    37.75,
    36.34139,
    36.56583,
    36.39111,
    35.85694,
    35.60472,
    35.68944000000001,
    35.44778,
    37.90222,
    36.69528,
    36.59444000000001,
    36.06528,
    35.66389,
    36.65139,
    35.39111,
    34.97694,
    35.18028,
    34.73028,
    35.00444,
    35.021390000000004,
    34.68639,
    34.69139000000001,
    34.68528,
    34.22611,
    35.50360999999999,
    35.47222,
    34.66167,
    34.396390000000004,
    34.18583,
    34.06583,
    34.34028,
    33.84167,
    33.55972,
    33.606390000000005,
    33.24944,
    32.74472,
    32.78972,
    33.23806,
    31.91111,
    31.56028,
    26.2125,
]

問題: 凸包

上の日本の都道府県の県庁所在地の経度・緯度のリストに対して, 凸包を構築せよ. 凸包上にある(すなわち日本の都道府県を包むような)点の番号を求めよ.

問題: Voronoi図

上の日本の都道府県の県庁所在地の経度・緯度のリストに対して, Voronoi図を描画せよ. この図の領域は何を意味する考察せよ.

問題: Euclid最小費用完全マッチング (難)

2次元平面 $[0,1]^2$ 上にランダムに分布した偶数個の点に対する最小費用の完全マッチング (すべての点の次数が $1$ の部分グラフ)を求める問題を考える. ただし枝の費用は点間のEuclid距離とする.このような問題をEuclid最小費用マッチング問題とよぶ. $K$-d木を用いて,各点に対して近い順に $10$個の点の間にだけ枝をひいた疎なグラフを作成し, それ上で最小費用マッチングを求める近似解法を考える. $100$点のEuclid最小費用マッチング問題をランダムに生成し, 疎なグラフ上で求めた近似マッチングと, すべての点間の枝を用いて最適解を求めた場合の費用と計算時間を測定せよ.

ヒント: 最小費用マッチング問題の求解には,NetworkXの関数max_weight_matchingを用いることができる. この関数は負の重みには対応していない. したがって,費用の合計を最小化するためには,大きな数(たとえば枝の費用の最大値)から枝の費用を減じた値を新たに枝の費用と 定義して最大化を行う必要がある.

問題:Euclid最短路 (難)

$(0,0)$ の地点から $(1,1)$ の地点にロボットが移動する経路を考える. ロボットは平面上のどの地点でも通過できるようにプログラムできるが, 平面 $[0,1]^2$ 上にランダムに配置された障害物を表す点の内側の領域は通過することができない. ランダムに $100$個の障害物を配置し,その外側を通過する最短距離を求めよ.

ヒント:$(0,0)$ と $(1,1)$ を加えた点に対する凸包を用いる.

確率・統計


確率・統計に関するサブモジュール stats

  • 確率分布の基礎
  • 共通のメソッド
  • 連続確率変数
  • 離散確率変数
  • データのあてはめ
  • 相関と回帰
  • 分布テスト
  • 信頼区間

例として平均 $100$,標準偏差 $10$ の正規分布 $N(100,10^2)$ の統計量を分布に付随するメソッドstatで表示してみる. 引数のmoments='mvsk'は,平均(mean),分散(variance),歪度(skewness),尖度(kurtosis)を意味する.

正規分布インスタンスを生成するnormは,位置パラメータlocで平均を,尺度パラメータscaleで標準偏差を与える.

密度関数 $f(x)$ をもつ連続な確率変数に対して, $X$ の期待値(平均)を $$ E[X]= \int_{-\infty}^{\infty} x f(x) dx $$ と定義する.

確率変数 $X$ と正数 $k$ に対する $k$次モーメント(moment)は,$m=E[X]$ としたとき, $E[ (X-m)^k ]$ で定義される.

分散(variance) $Var(X)$ は,2次モーメント $E[ (X-m)^2]$ と定義される. 分散は分布のばらつきを表す尺度である. $\sqrt{Var(X)}$ を標準偏差(standard deviation)とよび $\sigma$ と記す.

歪度(skewness)は $E[ (X-m)^3]/\sigma^3$,尖度(kurtosis)は $E[ (X-m)^4]/\sigma^4$ で定義され, それぞれ分布の非対称性ならびに尖りを表す尺度である.

正規分布の歪度は $0$ であり左右対称である.歪度が正の分布は右側の裾野が長くなり,逆に負の分布は左側の裾野が長くなる.

正規分布の尖度は $0$($3$とする場合もあるが,SciPyでは$0$と定義されている)であり, 尖度が正の分布は,正規分布と比べて山の頂上付近で尖って,裾野が長い形状であり, 尖度が負の分布は,頂上付近が丸く裾野が短い形状をもつ.

from scipy.stats import norm

print(norm(loc=100, scale=10).stats(moments="mvsk"))
(array(100.), array(100.), array(0.), array(0.))

分布の固定化(freezing)

  • パラメータを固定した分布インスタンスを生成する.
  • 生成した固定分布に対して,メソッドを用いて計算する.

たとえば,パラメータ(平均100,標準偏差100)を固定した正規分布インスタンスnを生成した後で,分布関数メソッドcdfを用いて 確率変数が $80$ 以下になる確率を求めるコードは,以下のようになる.

n = norm(loc=100, scale=100)
print(n.cdf(80))
0.42074029056089696

ガンマ分布(形状パラメータの使用例)

正のパラメータ $\lambda ,a$ て定義される分布

平均は $a/\lambda$,分散は $a/\lambda^2$

  • 形状 (shape)パラメータ a (分布によって名前は異なる;形状パラメータは必ず第1引数なので,名前を省略)
  • 尺度パラメータ scale($1/\lambda$ に対応)

例:ガンマ分布を生成して,平均と分散を計算

from scipy.stats import gamma

g = gamma(1, scale=10)  # a=1 (shape). 1/gamma=10 (scale)
print(g.stats(moments="mv"))  # mean =a/gamma, variance =1/gamma^2
(array(10.), array(100.))

共通メソッド


rvs: 確率変数にしたがう擬似乱数を生成するメソッド

  • 引数size: 生成する乱数の数(タプルを引数としたときは多次元配列)
from scipy.stats import norm

n = norm(loc=100, scale=10)
print("array=", n.rvs(size=5))
print("matrix=", n.rvs(size=(2, 3)))
array= [101.84651023 129.17878502 103.55253976 114.27602703  87.39188029]
matrix= [[ 95.0525987  110.6083437   76.65235316]
 [ 87.21302208  96.00269644  79.62090203]]

pdf: 連続確率変数の密度関数(probability density function)

任意の区間 $[a,b]$ に対して,確率変数 $X$ が区間内に入る確率 $P(X \in [a,b])$ が,

$$ P(X \in [a,b])=\int_{a}^{b} f(x) dx $$

と表されるとき,$f(x)$ を密度関数とよぶ. 密度関数で定義される連続な確率変数が $1$点をとる確率は $0$ である. SciPyにおいては,pdf(a)は $P(X \in [a,a+1])$ を返す.

pmf: 離散確率変数の確率関数(probability mass function)

確率関数 $p(x)$ は確率変数 $X$ が $x$ になる確率 $P(X=x)$ を表す.

cdf: 分布関数(cumulative distribution function)

(離散,連続)確率変数 $X$ の分布関数 $F(x)$ は $$ F(x) =P(X \leq x) $$ と定義される.

ppf: パーセント点関数(percent point function) (分布関数の逆関数)

  • 例:平均 100,標準偏差 10 の正規分布にしたがう需要をもつ商品に対して品切れ率を 5% にするための在庫量

    from scipy.stats import norm print( norm.ppf(0.95,loc=100,scale=10) )

    116.44853627

sf: 生存関数(survival function) ($1-F(x)$)

(離散,連続)確率変数 $X$ の生存関数は $$ P(X > x) $$ と定義される.

分布関数が与えた値以下の確率を返すのに対して,生存関数は与えた値より大きい確率を返す. (離散確率変数の場合には,「以上」でなく「より大きい」確率であることに注意する必要がある.)

  • 例: 平均 100,標準偏差 10 の正規分布にしたがう需要をもつ商品に対して,123個仕入れたときの品切れ確率

    from scipy.stats import norm print( norm.sf(123,loc=100,scale=10) )

    0.0107241100217

例題: 確率変数に対するメソッド

  1. 平均24,標準偏差8の正規分布にしたがう確率変数$X$が, $k=30$ 以下になる確率を求めよ.

  2. 同じ分布に対して,90%の確率で $k$以下になる $k$ を求めよ.

  3. 同じ分布に対して,$k=30$ 以上になる確率を求めよ.

  4. 同じ分布に対して,10%の確率で $k$以上になる $k$ を求めよ.

  5. 同じ分布に対して,20以上,28以下になる確率を求めよ.

  6. 同じ分布に対して, $[24-L,24+L]$ の間に入る確率が90%になるような $L$ を求めよ.

n = norm(loc=24, scale=8)
print(n.cdf(30))
print(n.ppf(0.9))
print(n.sf(30))  # もしくは 1- n.cdf(30)
print(n.ppf(1 - 0.1))  # 2と同じ
print(n.cdf(28) - n.cdf(20))  # 28以下の確率から20以下の確率を引く
print(24 - n.ppf(0.05))  # 確率が0.05(10%の半分)になるkを求め,平均24から引く
0.7733726476231317
34.2524125243568
0.2266273523768682
34.2524125243568
0.38292492254802624
13.158829015611783

expect(func) : 引数として与えられた関数funcの分布に対する期待値

標準正規分布の人口密度をもつ直線上の町に住む人たちが,町の中心にある町役場まで歩いて行くときの平均距離を求めてみよう. これは,標準正規分布に対して,引数を絶対値を表す関数として与えたときの期待値になる.

print(norm.expect(lambda x: abs(x)))
0.7978845608028651

例題: 損出関数

在庫理論では,品切れ量を計算するために確率変数 $x$ が $k$以上になるときの期待値を 用いる.これは損出関数(loss function)とよばれ,以下のように定義される. $$ G(k)= \int_{k}^{\infty} (x-k) f(x) dx $$

たとえば,$k=0.1$ のときの標準正規分布の損出関数は, 期待値の計算する際の下限lbを $0.1$ としたときの $x-k$ の期待値として以下のように計算される.

損出関数は,密度関数 $f(x)$ と生存関数 $1-F(x)$ を用いて,解析的に $G(k)= f(k) -k(1-F(k))$ と計算できる.

以下のコードから, 結果が一致していることが確認できる.

k = 0.1
print(norm.expect(lambda x: x - k, lb=k))
print(norm.pdf(k) - k * norm.sf(k))
0.35093533120471415
0.3509353312047147

問題: 在庫

平均 $100$,標準偏差 $10$ の正規分布にしたがう需要をもつ商品に対して,$120$個仕入れたとき, 品切れ費用が 1個あたり $100$円,在庫費用が 1個あたり $10$円としたときの期待費用を求めよ.

代表的な連続確率変数

uniform: 一様分布(uniform distribution)

$[\alpha,\beta]$ 上で定義された一様分布は密度関数 $1/(\beta-\alpha)$ をもつ. (SciPyでは離散な一様分布randintも使うことができる.)

位置パラメータlocは下限値 $\alpha$ を表し, 尺度パラメータscaleは上限値と下限値の差 $\beta-\alpha$ を表す. 平均は $(\beta-\alpha)/2$,分散は $(\beta-\alpha)^2 /12$ である.

from scipy.stats import uniform

u = uniform(loc=0, scale=10)
print("平均と分散 ", u.stats(moments="mv"))

x = np.linspace(-3, 15, 100)
y = u.pdf(x)
plt.plot(x, y)
平均と分散  (array(5.), array(8.33333333))

expon: 指数分布(exponential distribution)

指数分布とは,正のパラメータ $\lambda$ に対して,密度関数が $$ f(x) = \lambda e^{-\lambda x} \quad \forall x \geq 0 $$ で与えられる分布である.引数として指定する尺度パラメータscaleは $1/\lambda$ である. 分布関数は $$ F(x) = 1 - e^{-\lambda x} \quad \forall x \geq 0 $$ となる.平均は $1/\lambda$,分散は $1/\lambda^2$ である.

指数分布は無記憶性(memoryless property)をもつ唯一の連続分布 (離散分布では,後述する幾何分布が無記憶性をもつ.) として知られ, 稀な現象の時間間隔として用いられる

from scipy.stats import expon

ex = expon(scale=10)
print("平均と分散 ", ex.stats(moments="mv"))

x = np.linspace(-3, 100, 100)
y = ex.pdf(x)
plt.plot(x, y)
平均と分散  (array(10.), array(100.))

cauchy: Cauchy分布(Cauchy distribution)

Cauchy分布とは,分布の最頻値を表すパラメータ $x_0$, と広がりを表すパラメータ $\gamma$ に対して,密度関数が $$ \frac{1}{\pi\gamma\,\left[1 + \left(\frac{x-x_0}{\gamma}\right)^2\right]} $$ で与えられる分布である.

分布関数は, $$ \frac{1}{\pi} \arctan\left(\frac{x-x_0}{\gamma}\right)+\frac{1}{2} $$ である.Cauchy分布は平均をもたず,2次モーメントが無限大になる分布として知られている.

試しにSciPyで平均,分散,歪度,尖度を計算してみる.

from scipy.stats import cauchy

cau = cauchy(0, 1)
print(cau.stats(moments="mvsk"))

x = np.linspace(-10, 10, 100)
y = cau.pdf(x)
plt.plot(x, y)
(array(nan), array(nan), array(nan), array(nan))

平均,分散,歪度,尖度ともに,数値ではないことを表すnan(not a number)が表示された.

norm: 正規分布(normal distribution)

正規分布 $N(\mu,\sigma^2)$ は最も良く使われる連続分布であり, その密度関数は平均を表す位置パラメータ $\mu$ と標準偏差を表す尺度パラメータ $\sigma$ を用いて,

$$ f(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}} \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2} \right) $$

と記述される.

中心極限定理から(標準偏差をもつ)独立な分布の平均の分布は正規分布に近づくことがいえる. 特に,平均 $0$,標準偏差 $1$ の正規分布 $N(0,1)$ は標準正規分布とよばれ,その密度関数と分布関数は,それぞれ $$ \phi (x) = \frac{1}{\sqrt{2\pi}} \exp\!\left(-\frac{x^2}{2} \right) $$ と $$ \Phi (x) = \int_{-\infty}^x \phi (t) dt = \frac{1}{2} +\frac{1}{2} erf \left( \frac{x}{\sqrt{2}} \right) $$ で定義される.ここで $erf(x)$ は誤差関数であり, $$ erf(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt $$ である.

lognorm: 対数正規分布(log-normal distribution)

対数をとったときに正規分布になる分布である. 正規分布にしたがう確率変数が負の値も許すのに対して,対数正規分布にしたがう確率変数は正の値しかとらないという性質をもつ. また,独立同分布に従う確率変数の積は漸近的に対数正規分布にしたがうという,大数の法則の積バージョンの性質をもつ.

密度関数は位置パラメータ $\mu$ と形状パラメータ $\sigma$(正規分布と異なり尺度パラメータでないことを注意)を用いて, $$ f(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}} \exp\!\left(-\frac{(\ln x-\mu)^2}{2\sigma^2} \right) \quad \forall x >0 $$ と記述される.分布関数は標準正規分布の分布関数を $\Phi$ としたとき,$\Phi \left(\frac{\ln x - \mu}{\sigma} \right)$ となる.

平均は $e^{\mu+\frac{ \sigma^2 }{2} }$,分散は $e^{2\mu+\sigma^2}(e^{\sigma^2}-1)$ である.

from scipy.stats import lognorm

ln = lognorm(1, loc=5)
print(ln.stats(moments="mv"))

x = np.linspace(-1, 10, 100)
y = ln.pdf(x)
plt.plot(x, y)
(array(6.64872127), array(4.67077427))

halfnorm: 半正規分布(half-normal distribution)

$X$ が正規分布 $N(0,\sigma^2)$ のとき, $Y=|X|$ がしたがう分布を半正規分布とよぶ. 一般に平均が $0$ でない場合には,折り畳み正規分布 (folded-normal distribution)とよばれるが, SciPyでは両方とも halfnorm で生成できる.

密度関数は位置パラメータ $\mu$ と尺度パラメータ $\sigma$ を用いて, $$ f(x)= \frac{1}{\sqrt{2\pi\sigma^2}} \, e^{ -\frac{(x-\mu)^2}{2\sigma^2} } + \frac{1}{\sqrt{2\pi\sigma^2}} \, e^{ -\frac{(x+\mu)^2}{2\sigma^2} } $$ と記述される.分布関数は標準正規分布の誤差関数を $erf(x)$ としたとき, $$ \frac{1}{2}\left[ erf\left(\frac{x+\mu}{\sqrt{2\sigma^2}}\right) + erf\left(\frac{x-\mu}{\sqrt{2\sigma^2}}\right)\right] $$ となる.

平均は $$ \mu_Y = \sigma \sqrt{\frac{2}{\pi}} \,\, \exp\left(\frac{-\mu^2}{2\sigma^2}\right) + \mu \, erf\left(\frac{\mu}{\sqrt{2\sigma^2}}\right) $$ であり, 分散は $\mu^2 + \sigma^2 - \mu_Y^2$ である.

from scipy.stats import halfnorm

hn = halfnorm(loc=0, scale=np.sqrt(2) * 10)
print(ln.stats(moments="mv"))

x = np.linspace(-1, 100, 100)
y = hn.pdf(x)
plt.plot(x, y)
(array(6.64872127), array(4.67077427))

erlang: Erlang分布(Erlang distribution)

Erlang分布は,独立で同一の $n$ 個の指数分布の和の分布であり, パラメータ $\lambda>0$ と自然数 $n$ に対して,密度関数が $$ f(x)={\lambda^{n} x^{n-1} e^{-\lambda x} \over (n-1)!} \quad \forall x \geq 0 $$ で与えられる分布である.ここで $n$ は形状パラメータ,$1/\lambda$ は尺度パラメータである. 分布関数は $$ F(x)= 1 - \sum_{k=0}^{n-1} {(\lambda x)^{k} \over k!} e^{-\lambda x} \quad \forall x \geq 0 $$ となる.平均は $n/\lambda$,分散は $n/\lambda^2$ である.

Erlang分布は待ち行列理論に応用をもつ.

from scipy.stats import erlang

er = erlang(3, scale=10)
print(er.stats(moments="mv"))

x = np.linspace(-1, 100, 100)
y = er.pdf(x)
plt.plot(x, y)
(array(30.), array(300.))

gamma: ガンマ分布(gamma distribution)

ガンマ分布は,正のパラメータ $\lambda, a$ に対して,密度関数が $$ f(x) = \frac{ \lambda^a x^{a-1} e^{-\lambda x}}{\Gamma(a)} \quad \forall x \geq 0 $$ で与えられる分布である.ここで $\Gamma$ はガンマ関数である. $a$ が形状パラメータ,$\lambda$ が尺度パラメータである. 分布関数は $$ F(x) = \frac{\gamma(a, \lambda x)}{\Gamma(a)} \quad \forall x \geq 0 $$ となる.ここで $\gamma$ は不完全ガンマ関数である. 平均は $a/\lambda$,分散は $a/\lambda^2$ である.

ガンマ分布はErlang分布のパラメータ $n$ を実数値 $a$ に一般化した分布である.

from scipy.stats import gamma

ga = gamma(3.5, scale=10)
print(ga.stats(moments="mv"))

x = np.linspace(-1, 100, 100)
y = ga.pdf(x)
plt.plot(x, y)
(array(35.), array(350.))

logistic: ロジスティック分布(logistic distribution)

ロジスティック分布は,パラメータ $\mu, s$ に対して,密度関数が $$ f(x) = \frac{\exp(-\frac{x-\mu}{s})}{s(1+\exp(-\frac{x-\mu}{s}))^2} $$ で与えられる分布である.分布関数は $$ F(x) = \frac{1}{1 + e^{-(x - \mu)/s}} = \frac{1}{2}\left\{\tanh\left(\frac{x - \mu}{2 s}\right) + 1 \right\} $$ となる.平均は $\mu$,分散は $\pi^s \mu^2/3$ である.

ロジスティック分布は正規分布と同様に釣鐘型の密度関数と$S$字(シグモイド)型の分布関数をもつが, ロジスティック分布の方が裾野が長いという特徴をもつ. ロジスティック分布の尖度は $1.2$ であり,正規分布の $0$ と比べて大きいので,より裾野が長いことが分かる.

from scipy.stats import logistic, norm

lo = logistic(0, 1)
print(lo.stats(moments="mvsk"))

x = np.linspace(-10, 10, 100)
y = lo.pdf(x)
y2 = norm(0, 1).pdf(x)
plt.plot(x, y)
plt.plot(x, y2)
(array(0.), array(3.28986813), array(0.), array(1.2))

weibull_min: Weibull分布(Weibull distribution)

Weibull分布は正のパラメータ $c, \lambda$ に対して,密度関数が $$ f(x) = c \lambda \left(\lambda x \right)^{c-1} \exp(-(\lambda x)^{c}) \quad \forall x \geq 0 $$ で与えられる分布である.ここで $c$ は形状パラメータ,$1/\lambda$ は尺度パラメータである. 分布関数は $$ F(x)=1-\exp( -(\lambda x)^c) $$ であり,平均は $\Gamma (1+1/c)/\lambda$,分散は $\left(\Gamma (1+2/c)-\Gamma (1+1/c)\right)^2/\lambda^2$ である. ここで $\Gamma$ はガンマ関数である.

Weibull分布は故障現象や寿命を記述するために用いられる. パラメータ $c<1$ のとき故障率が時間とともに小さくなり,$c>1$ のとき時間とともに大きくなる. $c=1$ のとき故障率は一定であり,この場合には指数分布と一致する.

from scipy.stats import weibull_min

wei = gamma(1, scale=10)
print(wei.stats(moments="mv"))

x = np.linspace(-1, 100, 100)
y = wei.pdf(x)
plt.plot(x, y)
(array(10.), array(100.))

beta: ベータ分布(beta distribution)

ベータ分布は正のパラメータ $a, b$ に対して,密度関数が $$ f(x)= \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\, x^{a-1}(1-x)^{b-1} \quad \forall 0 \leq x \leq 1 $$ で与えられる分布である.ここで $\Gamma$ はガンマ関数であり,$a,b$ はともに形状パラメータである.

分布関数は $$ F(x)=I_x(a,b) $$ である.ここで $I_x$ は正則不完全ベータ関数である. 平均は $a/(a+b)$,分散は $ab/ \left( (a+b)^2 (a+b+1) \right)$ である.

ベータ分布は特定の区間の値をとる確率分布を表すときに有効であり,パラメータ $a,b$ を変えることによって様々な形状をとることができる. プロジェクト管理の手法として有名なPERT(Program Evaluation and Review Technique)においても確率的に変動する作業時間の近似として この分布を用いている.

from scipy.stats import beta

be = beta(2, 4)
print(be.stats(moments="mv"))

x = np.linspace(0, 2, 100)
y = be.pdf(x)
plt.plot(x, y)
(array(0.33333333), array(0.03174603))

t: $t$分布($t$ distribution)

$t$分布は正のパラメータ $\nu$ に対して,密度関数が $$ f(t) = \frac{\Gamma((\nu+1)/2)}{\sqrt{\nu\pi\,}\,\Gamma(\nu/2)} (1+t^2/\nu)^{-(\nu+1)/2} $$ で与えられる分布である.ここで $\Gamma$ はガンマ関数であり,$\nu$ は自由度(degree of freedom)とよばれる形状パラメータである. ちなみに $\nu$ の引数名はdegree of freedomの略でdfである.

分布関数は $$ F(t)=I_x(\nu/2,\nu/2) $$ である.ここで$I_x$ は正則不完全ベータ関数であり, $$ x = \frac{t+\sqrt{t^2+\nu}}{2\sqrt{t^2+\nu}} $$ とする. 平均は $0$,分散は $\nu>2$ のとき $\nu/(\nu-2)$,$1< \nu \le 2$ のとき無限大である.

いま $x_1,x_2,\ldots,x_n$ を平均 $\mu$ の独立な正規分布のサンプルとする. 標本平均と不偏分散を $$ \bar{x} = \frac{x_1+x_2+ \cdots+x_n}{n} $$ $$ s = \sqrt{ \frac{1}{n-1} \sum_{i=1}^n (x_i - \bar{x})^2 } $$ とする. このとき, $$ t = \frac{\bar{x} - \mu}{s/\sqrt{n}} $$ は,自由度 $\nu=n-1$ の $t$分布にしたがう. この性質を利用して,$t$分布は,標本値から母集団の平均値を統計的に推定する区間推定や,母集団の平均値の仮説検定に利用される.

平均との仮説検定はttest_1samp,独立な2つのサンプル間の仮説検定はttest_ind, 対応がある2つのサンプル間の仮説検定はttest_relで行うことができる.

ちなみに$t$分布は,$\nu=1$ のときCauchy分布になり,$\nu$ が無限大に近づくと正規分布に漸近する.

from scipy.stats import t

tdist = t(1)
print(lo.stats(moments="mvsk"))

x = np.linspace(-10, 10, 100)
y = tdist.pdf(x)
y2 = norm(0, 1).pdf(x)
plt.plot(x, y)
plt.plot(x, y2)
(array(0.), array(3.28986813), array(0.), array(1.2))

chi2: $\chi^2$分布($\chi^2$distribution)

$\chi^2$分布は正のパラメータ $k$ に対して,密度関数が $$ f(x)= \frac{1}{2^{\frac{k}{2}}\Gamma\left(\frac{k}{2}\right)}\; x^{\frac{k}{2}-1} e^{-\frac{x}{2}} \quad \forall x \geq 0 $$ となる分布である. ここで $\Gamma$ はガンマ関数であり, $k$ は自由度を表す形状パラメータである. $k$ の引数名はdfである.

分布関数は $$ F(x)= \frac{1}{\Gamma\left(\frac{k}{2}\right)}\;\gamma\left(\frac{k}{2},\,\frac{x}{2}\right) \quad \forall x \geq 0 $$ である.ここで $\gamma$ は不完全ガンマ関数である. 平均は $k$,分散は $2k$ である.

$x_1,x_2,\ldots,x_k$ を独立な標準正規分布のサンプルとするとき, $$ \sum_{i=1}^k (x_i)^2 $$ は自由度 $k$ の $\chi^2$分布にしたがう. この性質を利用して $\chi^2$分布は $\chi^2$検定chisquareやFriedman検定friedmanchisquareに用いられる.

from scipy.stats import chi2

ch = chi2(3)
print(ch.stats(moments="mv"))

x = np.linspace(0, 10, 100)
y = ch.pdf(x)
plt.plot(x, y)
(array(3.), array(6.))