SciPy (http://www.scipy.org/) は様々な科学技術計算のための実装を含むパッケージである．
ここでは，以下のサブパッケージについて解説する．
- optimize: 非線形最適化と方程式の根
- spatial: 計算幾何学
- stats: 確率・統計
- interpolate: 補間
- integrate: 積分
%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"))
複数変数用の minimize(f,x0)
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}
import plotly.graph_objs as go
from scipy.optimize import rosen
X, Y = np.mgrid[-2:2:0.05, -1:3:0.05]
Z = rosen([X, Y])
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=np.log(Z + 0.01)))
fig.update_layout(autosize=False, width=800, height=500)
fig.show()
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)
勾配の情報や2次導関数の情報 (Hessian)を使うと，誤差が減少することが見てとれる．
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.update_layout(autosize=False, width=800, height=500)
fig.show()
f = lambda x: (1.5 - x[0] + x[0]*x[1])**2 + (2.25 - x[0] + x[0]*x[1]**2)**2 + (2.625 - x[0] + x[0]*x[1]**3)**2
x0 = np.array([3., 4.])
res = minimize(f, x0, method="Nelder-Mead")
print(res.x)
res = minimize(f, x0, method="BFGS")
print(res.x)
f = (
lambda x: (x[0] - 10 * x[1]) ** 2
+ 5 * (x[2] - x[3]) ** 2
+ (x[1] - 2 * x[2]) ** 4
+ 10 * (x[0] - x[3]) ** 4
)
x0 = np.array([1.0, 1.0, 1.0, 1.0])
res = minimize(f, x0, method="Nelder-Mead")
print(res.x)
res = minimize(f, x0, method="BFGS")
print(res.x)
上の結果から，単体法（Nelder-Mead法）とBFGS（準Newton法）は，ほぼ同等の精度の解を算出することが分かる．
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.update_layout(autosize=False, width=800, height=500)
fig.show()
f = (
lambda x: -20 * np.exp(-0.2 * np.sqrt(0.5 * (x[0] ** 2 + x[1] ** 2)))
- np.exp(0.5 * (np.cos(2 * np.pi * x[0]) + np.cos(2 * np.pi * x[1])))
+ np.exp(1)
+ 20
)
x0 = np.array([1.0, 1.0])
res = minimize(f, x0, method="Nelder-Mead")
print(res.x)
res = minimize(f, x0, method="BFGS")
print(res.x)
上の結果から，多くの局所的最適解をもつ多峰性関数に対しては， 最適解の近くから探索をしないと， 最適解に到達しない場合があることが分かる．
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.update_layout(autosize=False, width=800, height=500)
fig.show()
制約付き最適化
半径 $\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)
問題（制約付き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$ 各家の水の消費量」の総和が最小になる場所に 井戸を掘ることにした． ただし, 各家の水の消費量は人数に比例するものとする．
どこを掘っても水が出るものとしたとき，どのようにして掘る場所を決めれば良いだろうか？
$(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]
import pandas as pd
import plotly.express as px
df = pd.DataFrame({"x": x, "y": y, "w": w})
fig = px.scatter(df, x="x", y="y", size="w", text=df.index)
fig.update_layout(autosize=False, width=600, height=500)
fig.show()
問題（ポートフォリオ最適化）
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}{lll} 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} $$問題（学区割り当て問題） （難）
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} $$数値は自分で適当に定めて解け．
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)
plt.plot(xdata, ydata, "ro")
plt.plot(xdata, f(xdata, param[0], param[1], param[2]));
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);
from scipy.optimize import root, brentq
print("Root finding ======\n ", root(f, x0=2.5))
print("Brentq=", brentq(f, a=1.0, b=3.0))
$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(1000, 2)
tree = KDTree(points)
dis, near = tree.query(x=[[0.1, 0.1], [0.5, 0.5]], k=10)
print(dis) # 各点から近い点への距離（近い順）
plt.plot(points[:, 0], points[:, 1], "b+")
plt.plot(points[near, 0], points[near, 1], "ro");
メソッド 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");
凸包クラス ConvexHull
凸包 (convex heull) とは，計算幾何学の基本となる概念であり，与えられたEuclid空間内の点を含む最小の凸多面体を指す．
凸包クラス ConvexHull のコンストラクタの引数は，以下のpoints配列である．
- 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)
Delaunay 3角形分割クラス Delaunay
Delaunay3角形分割 (Delaunay triangulation) とは， (2次元のときは）領域を点を通る3角形に分割したとき，どの3角形の頂点を通る円も，他の点を含まない3角形分割を指し． $K$-d木と同様に，近い点同士を見つけるのに便利であり，次に説明するVoronoi図と双対の関係にある．
Delaunay 3角形分割クラス Delaunay のコンストラクタの引数は，凸包と同じで点の配列である．
- 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(8, 2)
tri = Delaunay(points, furthest_site=False)
fig = delaunay_plot_2d(tri)
Voronoi図クラス Voronoi
Voronoi図 (Voronoi diagram)とは， 各点に近い空間で領域分けされた図であり，Delaunay3角形分割の双対グラフである．
与えられた点を母点，各母点に近い空間から成る領域を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)
tri = Delaunay(points)
fig = delaunay_plot_2d(tri)
Delaunay3角形分割の例題のVoronoi図を計算し，対応するDelaunay3角形分割と重ねて描画する．
np.random.seed(5)
points = np.random.rand(8, 2)
fig, ax = plt.subplots()
vor = Voronoi(points)
voronoi_plot_2d(vor, ax, show_vertices=False)
tri = Delaunay(points)
delaunay_plot_2d(tri, ax);
距離計算のためのサブパッケージdistance
距離計算のためのサブパッケージ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}$
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)
例題：Euclid最小木問題
点間の距離はEuclid距離の最小木問題を考える．ここで最小木とは，閉路を含まない連結グラフで，距離の合計が最小のものである． Euclid最小木問題の最適解に含まれる枝は，必ずVoronoi図の隣接する領域間の母点対の枝集合に含まれるという性質を利用すると，高速に解ける． Voronoi図の隣接する領域間の母点対は ridge_points を用いて列挙する． また，最小木を求めるアルゴリズムは，後述するNetworkXパッケージを用いる．
すべての枝をグラフに追加して最小木を求めたときと，Voronoi図を用いて必要な枝だけを追加して高速化した場合を比較する．
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))
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]
問題: Euclid最小費用完全マッチング （難）
2次元平面 $[0,1]^2$ 上にランダムに分布した偶数個の点に対する最小費用の完全マッチング （すべての点の次数が $1$ の部分グラフ）を求める問題を考える． ただし枝の費用は点間のEuclid距離とする．このような問題をEuclid最小費用マッチング問題とよぶ． $K$-d木を用いて，各点に対して近い順に $10$個の点の間にだけ枝をはった疎なグラフを作成し， それ上で最小費用マッチングを求める近似解法を考える． $100$点のEuclid最小費用マッチング問題をランダムに生成し， 疎なグラフ上で求めた近似マッチングと， すべての点間の枝を用いて最適解を求めた場合の費用と計算時間を測定せよ （ヒント： 最小費用マッチング問題の求解には，後述するNetworkXの関数max_weight_matchingを用いることができる． この関数は負の重みには対応していない． したがって，費用の合計を最小化するためには，大きな数（たとえば枝の費用の最大値）から枝の費用を減じた値を新たに枝の費用と 定義して最大化を行う必要がある）．
確率分布の基礎
例として平均 $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"))
通常は，パラメータを固定した分布インスタンスを生成し， 生成した固定分布に対して，メソッドを用いて計算する． これを分布の 固定化 (freezing)とよぶ，
たとえば，パラメータ（平均100，標準偏差100）を固定した正規分布インスタンス n を生成した後で，分布関数メソッド cdf を用いて 確率変数が $80$ 以下になる確率を求めるコードは，以下のようになる．
n = norm(loc=100, scale=100)
print(n.cdf(80))
分布によっては， （正規分布生成で用いた位置パラメータlocと尺度パラメータscaleの他に） 形状 (shape)パラメータが必要になる．
例として，ガンマ分布を生成する． この分布は， 正のパラメータ $\lambda ,a$ て定義され， 平均 $a/\lambda$，分散 $a/\lambda^2$ をもつ．
パラメータ $a$ が形状パラメータである． 形状パラメータは分布によって名前が異なるが， 形状パラメータは必ず第1引数なので，名前を省略して入力する． また， 尺度パラメータ scale は $1/\lambda$ に対応する．
$a=1, \gamma =1/10$ のガンマ分布を生成して，平均と分散を計算する．
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
from scipy.stats import norm
n = norm(loc=100, scale=10)
print("array=", n.rvs(size=5))
print("matrix=", n.rvs(size=(2, 3)))
from scipy.stats import norm
print(norm.ppf(0.95, loc=100, scale=10))
from scipy.stats import norm
print(norm.sf(123, loc=100, scale=10))
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から引く
print(norm.expect(lambda x: abs(x)))
例題： 損出関数
在庫理論では，品切れ量を計算するために確率変数 $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))
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)
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)
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)
平均，分散，歪度，尖度ともに，数値ではないことを表す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) $$ で定義される．ここで $\mathrm{erf}(x)$ は誤差関数であり， $$ \mathrm{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(4, 10, 100)
y = ln.pdf(x)
plt.plot(x, y);
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)
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)
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)
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)
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)
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)
t: $t$分布($t$ distribution)
$t$分布は正のパラメータ $\nu$ に対して，密度関数が $$ f(t) = \frac{\Gamma((\nu+1)/2)}{\sqrt{\nu\pi\,}\,\Gamma(\nu/2)} (1+ \frac{t^2}{\nu} )^{-(\nu+1)/2} $$ で与えられる分布である．ここで $\Gamma$ はガンマ関数であり，$\nu$ は自由度(degree of freedom)とよばれる形状パラメータである． ちなみに $\nu$ の引数名はdegree of freedomの略でdfである．
分布関数は $$ F(t)=I_x( \frac{\nu}{2}, \frac{\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)
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)
f: $F$分布 ($F$ distribution)
$F$分布は正のパラメータ $d_1, d_2$ に対して，密度関数が $$ f(x) = \frac{\sqrt{\frac{(d_1\,x)^{d_1}\,\,d_2^{d_2}} {(d_1\,x+d_2)^{d_1+d_2}}}} {x\,\mathrm{B}\!\left(\frac{d_1}{2},\frac{d_2}{2}\right)} \quad \forall x \geq 0 $$ で与えられる分布である．ここで $\mathrm{B}$ はベータ関数であり，$d_1, d_2$ は自由度を表す形状パラメータである． $d_1, d_2$ の引数名は，それぞれdfn,dfdである．
分布関数は $$ F(x)=I_{\frac{d_1 x}{d_1 x + d_2}} \left(\tfrac{d_1}{2}, \tfrac{d_2}{2} \right) \quad \forall x \geq 0 $$ である． 平均は $d_2>2$ のとき $d_2/(d_2-2)$，分散は $d_2>4$ のとき $$ \frac{2\,d_2^2\,(d_1+d_2-2)}{d_1 (d_2-2)^2 (d_2-4)} $$ である．
$U_1, U_2$ を独立な $\chi^2$分布にしたがう自由度 $d_1, d_2$ の確率変数としたとき， $$ \frac{U_1/d_1}{U_2/d_2} $$ は $F$分布にしたがうことが知られている． この性質を利用して$F$分布は $F$検定f_onewayに用いられる．
from scipy.stats import f
fdist = f(3, 3)
print(fdist.stats(moments="mv"))
x = np.linspace(0, 10, 100)
y = fdist.pdf(x)
plt.plot(x, y)
triang: 3角分布 (triangular distribution)
3角分布は，区間 $[\alpha,\beta]$ で $\gamma$ で最大になる密度関数をもった連続分布である． 密度関数は， $$ f(x)=\begin{cases} \frac{2(x-\alpha)}{(\beta-\alpha)(c-\alpha)} & \forall \alpha \le x < \gamma \\ \frac{2}{\beta-\alpha} & \forall x = \gamma \\ \frac{2(\beta-x)}{(\beta-\alpha)(\beta-\gamma)} & \forall c < x \le \beta \end{cases} $$ であり，分布関数は $$ F(x)=\begin{cases} \cfrac{(x-\alpha)^2}{(\beta-\alpha)(\gamma-\alpha)} &\forall \alpha \le x \le \gamma \\ 1-\cfrac{(\beta-x)^2}{(\beta-\alpha)(\beta-\gamma)} & \forall \gamma<x\le \beta \end{cases} $$ となる．
位置パラメータlocは下限値 $\alpha$ を表し， 尺度パラメータscaleは上限値と下限値の差 $\beta-\alpha$ を表し， 形状パラメータcは最頻値 $\gamma$ がloc+c*scaleになるように定められる． （つまり，形状パラメータcは $(\gamma-\alpha)/(\beta-\alpha)$ に設定する．）
平均は $$\frac{\alpha+b\beta\gamma}{3}$$ であり，分散は $$ \frac{\alpha^2+\beta^2+\gamma^2-\alpha \beta-\alpha \gamma-\beta \gamma}{12} $$ となる．
主に，最小値 $\alpha$ と最大値 $\beta$ と最頻値 $\gamma$ しか分かっていない場合に，簡易的に用いられる．
from scipy.stats import triang
tri = triang(0.2, loc=0, scale=10)
print(tri.stats(moments="mv"))
x = np.linspace(0, 10, 100)
y = tri.pdf(x)
plt.plot(x, y)
multivariate_normal: 多変量正規分布(multivariate normal distribution)
多変量正規分布は多次元の正規分布であり，平均を表す多次元ベクトル $\mu$ と共分散行列を表す尺度パラメータ $\Sigma$ を用いて， $$ f(x) = \frac{1}{\sqrt{(2 \pi)^k \det \Sigma}} \exp\left( -\frac{1}{2} (x - \mu)^T \Sigma^{-1} (x - \mu) \right) $$ と定義される．ここで $\det \Sigma$ は行列 $\Sigma$ の行列式を表す．
平均と共分散行列を固定化した分布を生成するには，引数meanとcovを用いる． 平均が $(0,0)$ で共分散行列が単位行列の2次元多変量正規分布を生成するには，以下のように記述する．
from scipy.stats import multivariate_normal
n = multivariate_normal([0.0, 0.0])
print(n.rvs(3))
print(n.pdf([0.0, 0.0]))
from scipy.stats import poisson
poi = poisson(4)
print(poi.stats(moments="mv"))
x = np.arange(0, 11)
y = poi.pmf(x)
plt.plot(x, y, marker="o")
geom: 幾何分布(geometric distribution)
幾何分布とは，パラメータ $0 < p <1$ に対して，確率関数が $$ P(X=k) = p (1-p)^{k-1} \quad \forall k=1,2,\ldots $$
で与えられる離散分布である．（この分布をファーストサクセス分布とよび，幾何分布を $p (1-p)^{k}$ と定義することもある．）
パラメータ $p$ は形状パラメータpとして入力する． 幾何分布は，独立なコイン投げ（表の出る確率を $p$ とする）を行ったときの表が出るまでに投げる回数を表す． 平均は $1/p$，分散は $(1-p)/p^2$ である．幾何分布は無記憶性 （$P(X >x+y | X>x)=P(X>y)$ がすべての $x,y=1,2,\ldots$ に対して成立すること） をもつ．
from scipy.stats import geom
ge = geom(0.5)
print(ge.stats(moments="mv"))
x = np.arange(0, 11)
y = ge.pmf(x)
plt.plot(x, y, marker="o")
binom: 2項分布(binomial distribution)
2項分布とは，パラメータ $0 < p <1$，自然数 $n$ に対して，確率関数が $$ P(X=k) = {n\choose k} p^k(1-p)^{n-k} \quad \forall k=0,1,\ldots,n $$ で与えられる離散分布である．ただし， $$ {n\choose k} = \frac{n!}{k! (n-k)!} $$ である．パラメータ $n,p$ は形状パラメータn,pとして入力する． 2項分布は，$n$回の独立なコイン投げ（表の出る確率を $p$ とする）を行ったときの表の出る数を表す． 平均は $np$，分散は $np(1-p)$ である．
from scipy.stats import binom
bi = binom(10, 0.3)
print(bi.stats(moments="mv"))
x = np.arange(0, 11)
y = bi.pmf(x)
plt.plot(x, y, marker="o")
nbinom: 負の2項分布(negative binomial distribution)
負の2項分布とは，パラメータ $p \leq 1$，正数 $n$ に対して，確率関数が $$ {k+n-1 \choose n-1} p^n (1-p)^k \quad \forall k=0,1,\cdots $$ で与えられる離散分布である．パラメータ $n,p$ は形状パラメータn,pとして入力する．
負の2項分布は，独立なコイン投げ（表の出る確率を $p$ とする）を行ったとき， $n$回表が出るまでに裏が出た回数を表す． 平均は $n(1-p)/p$，分散は $n(1-p)/p^2$ である．
from scipy.stats import nbinom
nbi = nbinom(4, 0.3)
print(nbi.stats(moments="mv"))
x = np.arange(0, 20)
y = nbi.pmf(x)
plt.plot(x, y, marker="o")
hypergeom: 超幾何分布(hypergeometric distribution)
超幾何分布とは，自然数 $M,n,N$ に対して，確率関数が $$ P(X=k) = \frac{\binom{n}{k}\binom{M-n}{N-k}}{\binom{M}{N}} \quad \forall \max (0, N - (M-n)) \leq k \leq \min (n, N) $$ で与えられる離散分布である． パラメータ $M, n, N$ は形状パラメータM, n, Nとして入力する． 超幾何分布は， $N$ 個の当たりをもつ $M$ 個入りのくじから， $n$ 個を非復元抽出したときに $k$ 個の当たりが含まれている確率を与える． 平均は $nN/M$，分散は $$ \frac{(M-n)n(M-N)N}{(M-1)M^2} $$ である．
from scipy.stats import hypergeom
hy = hypergeom(100, 10, 30)
print(hy.stats(moments="mv"))
x = np.arange(0, 20)
y = hy.pmf(x)
plt.plot(x, y, marker="o")
from scipy.stats import norm
rvs = norm(loc=100, scale=10).rvs(size=10)
print(norm.fit(rvs))
rvs = norm(loc=100, scale=10).rvs(size=10000)
print(norm.fit(rvs))
サンプル数が多い方が，良い推定値を返すことが確認できる．
平均や標準偏差を固定して推定することもできる．この場合には， 引数 floc, fscaleで位置パラメータと尺度パラメータの固定値を渡す．
例として，平均を $100$ に固定して推定してみる．
rvs = norm(loc=100, scale=10).rvs(size=100)
print(norm.fit(rvs, floc=100))
標準偏差の推定値が改善した．
ガンマ分布のように形状パラメータを必要とする場合には，引数名はf0, f1, ... と番号で指定する．
from scipy.stats import gamma
rvs = gamma(1, scale=10).rvs(size=10000)
print(gamma.fit(rvs))
形状パラメータ $a$ を $1$ に固定して予測する．
rvs = gamma(1, scale=10).rvs(size=10000)
print(gamma.fit(rvs, f0=1))
from scipy import stats
price = np.array(
[
7900,
11300,
13500,
6500,
5920,
6200,
6800,
1150,
8500,
7900,
8250,
5970,
20000,
7550,
6900,
10000,
1250,
]
)
t, p = stats.ttest_1samp(price, popmean=8000)
print(price.mean(), p)
speed = np.array([50, 29, 38, 43, 42, 40, 36, 37])
t, p = stats.ttest_1samp(speed, popmean=35)
print(speed.mean(), p / 2)
from scipy.stats import norm, probplot
from scipy.stats import uniform
u = uniform(loc=0, scale=10)
rvs = u.rvs(size=100)
xy_pair, corr = probplot(rvs, plot=plt)
print(corr)
正規分布に近いと基準線である赤線に近づく． 今度は，一様分布と近いか否かを引数 dist を uniform にして調べてみる．
xy_pair, corr = probplot(rvs, dist="uniform", plot=plt)
print(corr)
決定係数の平方根が $1$ に近づき，基準線と重なっているので，サンプルは正規分布より一様分布に近いと判断できる．
from scipy.stats import norm
import numpy as np
from scipy import stats
mu = 100.0
sigma = 10.0
normal = norm(mu, sigma)
population = normal.rvs(10000)
population
サンプルサイズ (sample size) $50$ のサンプルを， $1000$ 個（これをサンプル数と呼ぶ； number of samplesのこと．統計の用語はあまり良くない！）抽出し，その平均と標準偏差を計算する．
ただし，標準偏差は平均が未知であるので，不偏標準偏差を計算する必要がある． NumPyの標準偏差を計算するためのメソッド std では， 既定値は通常の（不偏ではない）標準偏差であり， 不偏にするためには引数 ddof（Delta Degrees of Freedomの略；自由度の差を表す）を $1$ に設定する必要がある．
$n$ 個のサンプル $x_1, x_2, \ldots, x_n$ があり、その平均値を $\bar{x}$ としたとき， 分散は以下のように計算される． $$ \sum_{i=1}^n (x_i-\bar{x})^2/n $$
一方，不偏分散は以下のように， 分母の $n$ から $1 (=$ ddof)を減じて補正したものである． $$ \sum_{i=1}^n (x_i-\bar{x})^2/(n-1) $$
平均は母集団の平均とほぼ同じで，標準偏差は母集団の標準偏差をサンプルサイズの平方根で割ったものとほぼ同じであることが確認できる（大数の法則）．
sample_size = 50
num_samples = 1000
sample_mean, sample_std = [], []
for i in range(num_samples):
sample = np.random.choice(population, sample_size)
sample_mean.append(sample.mean())
sample_mean = np.array(sample_mean)
print("標本平均と標本平均の標準偏差 = ", sample_mean.mean(), sample_mean.std(ddof=1))
print("母集団の標準偏差/sqrt(sample_size)= ", sigma / np.sqrt(sample_size))
１つのサンプルの標本平均 $\bar{x}$ と標本標準偏差 $s/\sqrt{n}$ （$n$ はサンプルサイズ）を用いて，90%信頼区間を計算する．
$100(1-\alpha)$%の信頼区間は，分布が正規分布で近似できる場合には，以下のようになる．
$$ \bar{x}-z_{\frac{\alpha}{2}}\sqrt{\frac{s^2}{n}} < \mu < \bar{x}+z_{\frac{\alpha}{2}}\sqrt{\frac{s^2}{n}} $$
ここで $\mu$ は母集団の平均（未知数）であり，上式が成立する確率が $100(1-\alpha)$%になるように信頼区間は設定される．
SciPyのstatsにある正規分布normのinterval(alpha, loc, scale)は， 平均loc，標準偏差scaleの正規分布の 100(1$-$alpha)%点を求めるメソッドであるので， サンプルの90%信頼区間は以下のように計算できる．
stats.norm.interval(0.9, sample.mean(), sample.std(ddof=1) / np.sqrt(sample_size))
サンプルをたくさん生成して信頼区間を計算し，母集団の平均 $\mu$ がその区間に含まれる割合を計算する．
おおよそ90%になっていることが確認できる．
count = 0
for i in range(num_samples):
sample = np.random.choice(population, sample_size)
lb, ub = stats.norm.interval(
0.9, sample.mean(), sample.std(ddof=1) / np.sqrt(sample_size)
)
if mu > lb and mu < ub:
count += 1
print(count / num_samples)
今度は，サンプルサイズを $10$ に設定してみる．一般に，サンプルサイズ $n$ が $30$ 未満の場合には，正規分布で近似するのではなく，$t$ 分布を用いることが推奨される．
$100(1-\alpha)$%の信頼区間は，分布が$t$ 分布で近似できる場合には，以下のようになる．
$$ \bar{x}-t_{\frac{\alpha}{2}}(n-1)\sqrt{\frac{s^2}{n}} < \mu < \bar{x}+t_{\frac{\alpha}{2}}(n-1)\sqrt{\frac{s^2}{n}} $$
まずは正規分布と仮定して，信頼区間を計算して，母集団平均 $\mu$ が区間に含まれている割合を計算してみる． 指定した90%ではなく，それより小さい値になっていることが確認できる．
sample_size = 10
num_samples = 1000
sample_mean, sample_std = [], []
for i in range(num_samples):
sample = np.random.choice(population, sample_size)
sample_mean.append(sample.mean())
stats.norm.interval(0.9, sample.mean(), sample.std(ddof=1) / np.sqrt(sample_size))
count = 0
for i in range(num_samples):
sample = np.random.choice(population, sample_size)
lb, ub = stats.norm.interval(
0.9, sample.mean(), sample.std(ddof=1) / np.sqrt(sample_size)
)
if mu > lb and mu < ub:
count += 1
print(count / num_samples)
$t$ 分布を用いて信頼区間を計算する． t.intervalの第２引数は自由度であり， サンプルサイズ $-1$ と設定する．
正規分布より広い区間となっており，母集団平均もほぼ90%の割合で区間内に入っていることが確認できる．
stats.t.interval(
0.9, len(sample) - 1, sample.mean(), sample.std(ddof=1) / np.sqrt(sample_size)
)
count = 0
for i in range(num_samples):
sample = np.random.choice(population, sample_size)
lb, ub = stats.t.interval(
0.9, len(sample) - 1, sample.mean(), sample.std(ddof=1) / np.sqrt(sample_size)
)
if mu > lb and mu < ub:
count += 1
print(count / num_samples)
例として1次元の関数 $f(x)=\sin x/(1+x^2)$ の補間を考える．
関数 $f(x)$ と$[-4,4]$ の整数に対応する $9$個のデータ点 $x_i,y_i=f(x_i)$ を生成しておく．
from scipy import interpolate
x = np.arange(-4, 5)
y = np.sin(x) / (1 + x**2)
plt.plot(x, y, "o");
Lagrange多項式
$k+1$ 個のデータ点 $(x_0,y_0), (x_1,y_1),\ldots,(x_k,y_k)$ が $x_i$ の昇順に並んでいるものとする
補間は，すべての点を通る多項式を1つ決めるのか， 区分 $[x_i, x_{i+1}]$ ごとに異なる多項式を用いるのかによって大きく2つに分けられる．
前者の方法の1つとして，Lagrange多項式(Lagrange polynomial)がある．
Lagrange基底多項式 $$ \ell_j(x) = \frac{(x-x_0)}{(x_j-x_0)} \cdots \frac{(x-x_{j-1})}{(x_j-x_{j-1})} \frac{(x-x_{j+1})}{(x_j-x_{j+1})} \cdots \frac{(x-x_k)}{(x_j-x_k)} $$ を用いて $$ L(x) = \sum_{j=0}^{k} y_j \ell_j(x) $$ と定義される．$\ell_j(x_j)=1$ であり，かつ $\ell_j(x_i)=0$ $(i \neq j)$ であるので， $L(x)$ はすべてのデータ点を通過することが分かる．
interpolateパッケージにあるlagrange(x,y)は1次元データに対するLagrange多項式を求めるための関数である．
これは，同じ長さをもつ1次元配列 $x$ と $y$ に対して， 関数 $y=f(x)$ を近似するLagrange多項式を返す． 返値はNumPyの1次元の多項式オブジェクトである．
f = interpolate.lagrange(x, y)
print(f)
xnew = np.arange(-4, 4, 0.01)
ynew = f(xnew)
plt.plot(x, y, "o", xnew, ynew, "-")
線形補間
区分ごとに異なる多項式を用いる場合には，自由度が増えるので様々な方法が考えられる．
interpolateパッケージにあるinterp1d(x,y)は1次元データに対する補間を行う関数である． 同じ長さをもつ1次元配列 $x$ と $y$ に対して， 関数 $y=f(x)$ を近似する補間を行う． 返値は，関数のように与えられた新しい点xに対して値yを返すオブジェクトである． $x,y$ 以外の主な引数は以下の通り．
kind: 補間の種類を表す文字列であり， "linear"（線形補間；既定値）， "nearest"（最近点補間）， "zero"（$0$次補間）， "slinear"（スプライン線形補間）， "quadratic"（スプライン2次補間）， "cubic"（スプライン3次補間）から選択する．
copy: 与えられたデータ $x,y$ のコピーを作成してから処理するか否かを表す論理値である．既定値はコピーすることを表すTrueである．
fill_value: 領域外の点を与えたときに割り当てられる値である．既定値はNaNである．
assume_sorted: 点を表すデータ $x$ が昇順に並んでいるか否かを表す論理値である．既定値は昇順に並んでいるデータを仮定するTrueである．
最も簡単な方法は区分ごとに異なる線分で繋いだ（区分的線形）関数である． これは線形補間(linear interpolation)とよばれる．
f = interpolate.interp1d(x, y, kind="linear")
xnew = np.arange(-4, 4, 0.01)
ynew = f(xnew)
plt.plot(x, y, "o", xnew, ynew, "-");
f = interpolate.interp1d(x, y, kind="nearest")
xnew = np.arange(-4, 4, 0.01)
ynew = f(xnew)
plt.plot(x, y, "o", xnew, ynew, "-");
f = interpolate.interp1d(x, y, kind="zero")
xnew = np.arange(-4, 4, 0.01)
ynew = f(xnew)
plt.plot(x, y, "o", xnew, ynew, "-");
f = interpolate.interp1d(x, y, kind="slinear")
xnew = np.arange(-4, 4, 0.01)
ynew = f(xnew)
plt.plot(x, y, "o", xnew, ynew, "-");
f = interpolate.interp1d(x, y, kind="quadratic")
xnew = np.arange(-4, 4, 0.01)
ynew = f(xnew)
plt.plot(x, y, "o", xnew, ynew, "-");
f = interpolate.interp1d(x, y, kind="cubic")
xnew = np.arange(-4, 4, 0.01)
ynew = f(xnew)
plt.plot(x, y, "o", xnew, ynew, "-");
問題： テーマパークの待ち時間 （欠損あり）
再び，アトラクションの待ち時間の調査を行った． 待ち時間を $30$ 分おきに測定した結果は， 開園時刻を $0$ としたPythonのリストとして，以下のように与えられている． 開園から閉園までの任意の時刻における待ち時間をLagrange補間と3次スプライン補間を用いて推定せよ． ただし，リスト内の NaN はNot a Number（数でないもの）を意味し，測定機械のミスでデータが取れなかったことを表す．
[5,55,80,60,70,60,50,60,70,55,40,55,55,NaN,NaN,80,80,80,70,80,80,70,60,45,30,25,20,15]
from scipy import integrate
f = lambda x: np.pi * x**2
print(integrate.quad(f, 0.0, 10.0)) # 積分
print(np.pi / 3.0 * 10**3) # 理論値
2重積分関数 dblquad
dblquadはFortran言語で書かれたQUADPACKを用いて2重定積分を求める．
引数は以下の通り．
- func: 少なくとも2つの変数（たとえば $x,y$）をもつ関数 $f(x,y)$ である．
- a: 積分範囲 $[a,b]$ の下限 $a$ を表す．$-\infty$ にしたい場合には，NumPyの-infを用いる．
- b: 積分範囲 $[a,b]$ の上限 $b$ を表す． $b> a$ である必要がある．$+\infty$ にしたい場合には， NumPyのinfを用いる．
- gfun: 第2引数である変数 $y$ の積分範囲の下限を規定する関数 $g(x)$ である．
- hfun: 第2引数である変数 $y$ の積分範囲の上限を規定する関数 $h(x)$ である．
返値は積分値と誤差のタプルである．
例題：正方形の都市
面積 $1$ の正方形の都市に一様に住んでいる人たちが，正方形の中心にあるスーパーに行くときの距離の期待値を求めてみよう． 正方形 $[0,1]^2$ 内の座標 $x,y$ に住んでいる人の中心への距離は $$ \sqrt{ \left(x- \frac{1}{2} \right)^2 + \left(y- \frac{1}{2} \right)^2 } $$ と計算できるので，これを区間 $[0,1]^2$ で積分することによって，期待値が計算できる． $x$ の積分範囲は $[0,1]$ で良いが，$y$ の積分範囲は $x$ の関数として与える必要があるので，ダミーの関数g,hを作成して 引数として渡していることに注意されたい．
理論値 $( \sqrt{2} + \log (1+\sqrt{2}) )/6 $ の導出は結構面倒である（拙著「はじめての確率論」（近代科学社） p.33 参照）．
f = lambda x, y: np.sqrt((x - 0.5) ** 2 + (y - 0.5) ** 2)
g = lambda x: 0.0
h = lambda x: 1.0
print(integrate.dblquad(f, 0, 1, g, h)) # 2重積分
print((np.sqrt(2) + np.log(1 + np.sqrt(2))) / 6.0) # 理論値