カテゴリー: 教養

MNISTから手書き文字認識用のデータセットをロードするクラスを作ってみた。 ロードしたデータセットをpickleでシリアライズ、デシリアライズする機能付き。 後々改造する予定でここに貼ったのはメンテしない予定。 たったこれだけ書くのに40分も要してしまった...。 Pythonの冗長感が半端ないけども慣れるしかない。 マジックコードだらけだけども、MNISTの手書き文字認識用データ取得専用だから仕方ない。 get_image()関数により訓練データ、訓練ラベル、テストデータ、テストラベルを取れる。 データは1次元になって入る。つまり1ファイルごとに28*28のサイズがある。 import os import urllib.request import numpy as np import gzip import pickle class MnistLoader: def __init__(self, mnistdir): self.url_base = \'http://yann.lecun.com/exdb/mnist/\' self.dataset_dir = mnistdir self.save_path = self.dataset_dir + \'mnist.pkl\' self.key_file = { \'train_img\':\'train-images-idx3-ubyte.gz\', \'train_label\':\'train-labels-idx1-ubyte.gz\', \'test_img\':\'t10k-images-idx3-ubyte.gz\', \'test_label\':\'t10k-labels-idx1-ubyte.gz\' } self.dataset = {} if os.path.exists(self.save_path): with open(self.save_path,\'rb\') as f: self.dataset = pickle.load(f) else: self.__load_mnist() self.dataset[\'train_img\'] = self.__load_img(self.key_file[\'train_img\']) self.dataset[\'train_label\'] = self.__load_label(self.key_file[\'train_label\']) self.dataset[\'test_img\'] = self.__load_img(self.key_file[\'test_img\']) self.dataset[\'test_label\'] = self.__load_label(self.key_file[\'test_label\']) with open(self.save_path, \'wb\') as f: pickle.dump(dataset, f, -1) def __load_mnist(self): \'\'\' load mnist data and store to file \'\'\' for v in self.key_file.values(): file_path = self.dataset_dir + v urllib.request.urlretrieve(url_base + v, self.dataset_dir + v) def __load_img(self, file_name): file_path = self.dataset_dir + file_name with gzip.open(file_path, \'rb\') as f: data = np.frombuffer(f.read(), np.uint8, offset=16) data = data.reshape(-1, 784) return data def __load_label(self,file_name): file_path = dataset_dir + \'/\' + file_name with gzip.open(file_path, \'rb\') as f: labels = np.frombuffer(f.read(), np.uint8, offset=8) return labels def get_image(self): return (self.dataset[\'train_img\'] , self.dataset[\'train_label\']), (self.dataset[\'test_img\'], self.dataset[\'test_label\']) mnist_loader = MnistLoader(mnistdir=\'/Users/ikuty/Documents/mnist/\') PILを使って訓練データの1枚目を表示するテスト。 import sys, os sys.path.append(os.pardir) import numpy as np from PIL import Image def img_show(img): pil_img = Image.fromarray(np.uint8(img)) pil_img.show() (x_train, t_train), (x_test, t_test) = mnist_loader.get_image() img = x_train[0] label = t_train[0] bimg = img.reshape(28, 28) img_show(bimg) こんなのが出る。(28x28しかなくて小さすぎなので256x256に引き伸ばして表示。)

[mathjax] 出力層の総和を1にするように調整できれば、出力層を確率としてとらえることができるようになる。 入力層に画像を放り込み、出力層でラベルに属する確率を出せば、画像にラベルをつける分類器になる。 入力が(n)ノードあるとして、全てのノードを入力を受けて出力する。 出力側のノードが(m)ノードあるとして、和が(1)になるようにどう操作するか。 分母は入力側の全ノードの和、分子は対応する1ノードだけ。そうやって1にする。 指数関数に入力を食わせて和を取る。指数関数は単調増加関数だから、 入力の大小と出力の大小は一致する。定義上は以下の通り。 begin{eqnarray} y_k &=& frac{e^{a_k}}{sum_{i=1}^n e^{a_i}} end{eqnarray} import numpy as np def softmax1(x): return np.exp(x) / np.sum(np.exp(x)) a = np.array([1010,1000,990]) r1 = softmax1(a) r1 # array([nan, nan, nan]) ただ、この通りにするとオーバーフローを起こすので、 式変形してオーバーフローを回避する。(C)は入力ノードの最大値。 begin{eqnarray} y_k &=& frac{e^{a_k}}{sum_{i=1}^n e^{a_i}} \\ &=& frac{C e^{a_k}}{Csum_{i=1}^n e^{a_i}} \\ &=& frac{e^{a_k+log{C}}}{sum_{i=1}^{n} e^{a_i + log{C}}} \\ &=& frac{e^{a_k+C\'}}{sum_{i=1}^{n} e^{a_i + C\'}} end{eqnarray} import numpy as np def softmax2(x): c = np.max(x) return np.exp(x-c) / np.sum(np.exp(x-c)) a = np.array([1010,1000,990]) r2 = softmax2(a) r2 # array([9.99954600e-01, 4.53978686e-05, 2.06106005e-09]) SoftMax後の和は1 前述の通り、SoftMax後の和は1。 z = np.sum(r2) z # 1.0

[mathjax] 順伝搬(forward)はNumpyの積和計算を使って超絶簡単に記述できる。 行列の積和計算 何はともあれNumpyで行列の積を計算する方法について。 この仕組みがあるから、forなどの制御構造を使わずに行列の積を実現できる。 begin{eqnarray} A &=& begin{pmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} end{pmatrix} \\ B &=& begin{pmatrix} b_{11} & b_{21} \\ b_{12} & b_{22} end{pmatrix} \\ A cdot B &=& begin{pmatrix} a_{11} b_{11} + a_{21} b_{12} && a_{11} b_{21} + a_{21} b_{22} \\ a_{12} b_{11} + a_{22} b_{12} && a_{12} b_{21} + a_{22} b_{22} end{pmatrix} end{eqnarray} 積はNumpyのdot関数で一発。 左行列と右行列の\"対応する行\"は一致している必要がある。 つまり(N_1 times M_1 cdot N_2 times M_2)という積計算において(M_1 = N_2)である必要がある。 a11,a21,a12,a22 = (1,2,3,4) b11,b21,b12,b22 = (1,2,3,4) a = np.array([[a11,a21],[a12,a22]]) b = np.array([[b11,b21],[b12,b22]]) np.dot(a,b) # array([[ 7, 10], # [15, 22]]) 和はNumpyの+オペレータで一発。 begin{eqnarray} A &=& begin{pmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} end{pmatrix} C &=& begin{pmatrix} c_{1} \\ c_{2} end{pmatrix} end{eqnarray} begin{eqnarray} A+C &=& begin{pmatrix} a_{11} + c_1 & a_{21} + c_1 \\ a_{12} + c_1 & a_{22} + c_2 end{pmatrix} end{eqnarray} c1, c2 = (5,5) c = np.array([[c1],[c2]]) a + c # array([[6, 7], # [8, 9]]) 各層における信号伝達の実装 Numpyで積和が簡単にかける。 すなわち、各層における伝達が積和で表されるとすると、 入力層から各層における信号伝達の実装を簡単にかける。 左側の層が(X=begin{pmatrix}x_1 x_2end{pmatrix})、重みが(W=begin{pmatrix}w_{11} & w_{21} \\ w_{12} & w_{22} end{pmatrix})、 右側の層が(A=begin{pmatrix} a_{11} & a_{21} end{pmatrix})とする。 すると、右層と左層の関係は以下の通り。 begin{eqnarray} A &=& XW \\ begin{pmatrix} a_{11} & a_{21} end{pmatrix} &=& begin{pmatrix}x_1 x_2end{pmatrix} begin{pmatrix}w_{11} & w_{21} \\ w_{12} & w_{22} end{pmatrix} \\ &=& begin{pmatrix} w_{11}x_1 + w_{21} x_2 \\ w_{12}x_1 + w_{22}x_2end{pmatrix} end{eqnarray} 左層にバイアス(B=begin{pmatrix}b_1 \\ b_1 end{pmatrix})があった場合、線形和で表現できる。 行と列がどちらの方向に対応するかに注意。 1つの層に(N)個のノードが存在することを(Ntimes 1)の行列で表現する。 begin{eqnarray} A &=& XW + B \\ begin{pmatrix} a_{11} & a_{21} end{pmatrix} &=& begin{pmatrix}x_1 x_2end{pmatrix} begin{pmatrix}w_{11} & w_{21} \\ w_{12} & w_{22} end{pmatrix} + begin{pmatrix}b_1 \\b_1 end{pmatrix} \\ &=& begin{pmatrix} w_{11}x_1 + w_{21} x_2 + b_1 \\ w_{12}x_1 + w_{22}x_2 + b_1end{pmatrix} end{eqnarray} x1,x2 = (0.4,0.9) w11,w21,w12,w22 = (0.8,0.5,0.1,0.3) b1 = 0.1 x = np.array([x1,x2]) w = np.array([[w11,w21],[w12,w22]]) b = np.array([b1,b1]) # これで一発 a = np.dot(x,w) + b # array([0.51, 0.57]) # 活性化関数としてsigmoid関数をかます z = sigmoid(a) # array([0.62480647, 0.63876318]) 順伝播 入力層から出力層まで上記のような計算を繰り返していく。これを\"順伝播\"とか言う。 順伝\"搬\"ではなく、順伝\"播\"。ネットワークが多層になったとしても全く同様。 以下みたいにNumpyの積和計算を繰り返すだけで出来る。 def init_network(): network = {} network[\'W1\'] = np.array([[0.1,0.3,0.5],[0.2,0.4,0.6]]) network[\'b1\'] = np.array([0.1,0.2,0.3]) network[\'W2\'] = np.array([[0.1,0.4],[0.2,0.5],[0.3,0.6]]) network[\'b2\'] = np.array([0.1,0.2]) network[\'W3\'] = np.array([[0.1,0.3],[0.2,0.4]]) network[\'b3\'] = np.array([0.1,0.2]) return network def forward(network, x): W1, W2, W3 = network[\'W1\'], network[\'W2\'], network[\'W3\'] b1,b2,b3 = network[\'b1\'], network[\'b2\'], network[\'b3\'] a1 = np.dot(x, W1) + b1 z1 = sigmoid(a1) a2 = np.dot(z1, W2) + b2 z2 = sigmoid(a2) a3 = np.dot(z2, W3) + b3 y = a3 return y network = init_network() x = np.array([1.0,0.5]) y = forward(network,x) # array([0.31682708, 0.69627909])

[mathjax] 深層学習入門。Python,Numpyにも少し慣れてきたので、 Numpyだけで伝搬,逆伝搬を計算することで深層学習に慣れていく。 単に自分の理解のためだけの記事なので、誤りがあっても気にしない。 活性化関数 まず活性化関数。activation function。 3つ(x_1,x_2,x_3)の入力があったとする。それぞれに重み(w_1,w_2,w_3)がかかるとする。 ノードは入力に重みをかけた和、つまり(w_1x_1 + w_2x_2 + w_3x_3)を受けるものとする。 ノードは受けた値の大きさに応じて出力を返す機能を持つ。 受けた値と出力の関係（まさに関数）を活性化関数と言って、 ステップ関数,Sigmoid関数,ReLU,Softmax関数などいくつか種類がある。 Step関数,Sigmoid関数,ReLUをNumpyだけで実装してみる。 import numpy as np import matplotlib.pyplot as plt # sigmoid function def sigmoid(x): return 1 / (1+np.exp(-x)) # step function def step(x): return np.array(x > 0, dtype=int) # ReLU def relu(x): return np.maximum(0,x) / 5 x = np.arange(-5.0,5.0,0.1) y = sigmoid(x) x = np.arange(-5.0,5.0,0.1) y1 = sigmoid(x) y2 = step(x) y3 = relu(x) plt.plot(x,y1) plt.plot(x,y2) plt.plot(x,y3) plt.ylim(-0.1,1.1) plt.show() どんなに入力信号が大きくても、出力を0から1の間に押し込める。 入力が大きければ出力が大きいという意図はあるものの、 入力と出力の関係が非線形になっているものが多い。 活性化関数が線形関数だと、ネットワークを重ねていったとしても、 ネットワークの最初の入り口と最後の出口を見たとして、1つの線形関数を通しただけ... ということになる。つまり加算と定数倍は何回実行したとしても、 一つの定数倍、加算の計算にまとめることができる。 これだと、層を重ねる意味がなくなってしまう。 活性化関数が非線形だと、重ねた活性化関数をまとめることはできず、 複雑な入出力を表現できるため、活性化関数として非線形関数を使用する。

[mathjax] マルコフ連鎖モンテカルロ法。2変量正規分布からGibbs Samplingする方法を考えてみました。 式を流してもよくわからないので、行間ゼロで理解できるまで細切れにして書いてみます。 Gibbs Sampling 無茶苦茶一般的に書かれている書籍だと以下みたいになっている。 ステップごとに(theta=(theta_1,theta_2,cdots,theta_n))の各要素を、 その要素以外の要素の条件付き確率でランダムに発生させる。 begin{eqnarray} theta^0 &=& (theta_1^0,theta_2^0,cdots,theta_n^0) \\ theta_i^{t+1} &sim& p(theta_i|theta_1^{t+1},cdots,theta_{i-1}^{t+1},theta_{i+1}^t,cdots,theta_n^t) end{eqnarray} 2次元空間であれば、ある点を決める際に、 片方の変数(theta_1)を固定して条件付き確率(P(theta_2|theta_1))を最大にするパラメタ(hat{theta_2})を決める。 その時点で((theta_1,hat{theta_2}))が決まる。 次は変数(hat{theta_2})を固定して条件付き確率(P(theta_1|hat{theta_2}))を最大にするパラメタ(hat{theta_1})を決める。 2次元正規分布であれば、片方の確率変数を固定したときの条件付き確率が1変数の正規分布になるから、 これがすごくやりやすい。 Gibbs Samplingは、このように条件付き確率を計算できる必要がある。 2変量正規分布の条件付き確率 2変量正規分布の片方の確率変数を固定すると1変数の正規分布が出てきます。 山の輪切りにして出てきた正規分布の母数を調べたいのですが、導出が無茶苦茶大変そうです。 出てきた正規分布の母数について式変形すると出てくるようです。こちらを参考にさせて頂きました。 この関係式を利用してひたすら式変形していきます。 begin{eqnarray} f(x|y) &=& frac{f(x,y)}{f(x)} end{eqnarray} 今、確率変数(y=y\')を固定し確率変数(x)の正規分布を考えます。 2変量正規分布の分散共分散行列が(sum)であるとします。 begin{eqnarray} sum &=& begin{pmatrix} sigma_{xx} & sigma_{xy} \\ sigma_{yx} & sigma_{yy} end{pmatrix} end{eqnarray} 出てくる正規分布の平均は以下となります。 確率変数が(x)なのに固定した(y)が出てくるのがポイント。 begin{eqnarray} mu\' = bar{x} + frac{sigma_{xy}}{sigma_{yy}} (y\'-bar{y}) end{eqnarray} また、出てくる正規分布の分散は以下となります。 begin{eqnarray} sigma\'^2 &=& sigma_{xx} - frac{sigma_{xy}}{sigma_{yy}} sigma_{yx} end{eqnarray} Python実装例 実際にPythonコードにしてみた図です。2変量正規分布の条件付き確率さえわかってしまえば あとはコードに落とすだけです。 import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats n_dim = 2 def gibbs_sampling(mu, sigma, sample_size): samples = [] start = [0, 0] samples.append(start) search_dim = 0 for i in range(sample_size): search_dim = 0 if search_dim == n_dim-1 else search_dim + 1 prev_sample = samples[-1][:] s11 = sigma[search_dim][search_dim] s12 = sigma[search_dim][search_dim - 1] s21 = sigma[search_dim -1 ][search_dim] s22 = sigma[search_dim - 1][search_dim - 1] mu_x = mu[search_dim] mu_y = mu[search_dim-1] _y = prev_sample[search_dim - 1] new_mean = mu_x + s12/float(s22) * (_y - mu_y) new_sigma = s11 - s12/float(s22) * s21 sample_x = np.random.normal(loc=new_mean, scale=np.power(new_sigma, .5), size=1) prev_sample[search_dim] = sample_x[0] samples.append(prev_sample) return np.array(samples) nu = np.ones(2) covariance = np.array([[0.5, 0.5], [0.5, 3]]) samples = gibbs_sampling(nu, covariance, 1000) fig, ax1 = plt.subplots(figsize=(6, 6)) ax1.scatter(sample[:, 0], sample[:, 1], marker=\"o\", facecolor=\"none\", alpha=1., s=30., edgecolor=\"C0\", label=\"Samples\" ) 実行結果

[mathjax] 最尤推定とMAP推定とベイズの定理は繋がっていたので、 記憶が定かなうちに思いの丈を書き出してみるテスト。俯瞰してみると面白い。 あるデータ達(x)が観測されていて、それらは未知のパラメータを持つ確率分布から発生している。 観測されたデータ達(x)を使って、それらのデータを発生させたモデルのパラメータを推定したい。 確率密度関数の中に2つの変数があって。 片方を定数、片方を確率変数として扱うことで2通りの見方ができる。 例えば(n)回のコイントスで(k)回表が出る確率が(theta)だとしてベルヌイ分布の確率密度関数は (k)と(theta)のどちらが確率変数だとしても意味がある。 begin{eqnarray} f(k;theta)=theta^k (1-theta)^{n-k} end{eqnarray} 表が出る確率(theta)が定数だと思って、確率変数(x)の確率密度関数と思う。 単なる確率変数(x)の確率密度関数の中に(theta)という定数がある。尤度。 尤度は確率変数(x)の確率密度関数!。 begin{eqnarray} p(X=x|theta) end{eqnarray} 尤度(p(x|theta))を最大にする(theta)を推定するのが最尤推定。 begin{eqnarray} newcommand{argmax}{mathop{rm arg~max}limits} hat{theta} = argmax_{theta} p(X=x|theta) end{eqnarray} 事後確率と事前確率には関係があって以下のようになる。ベイズの定理。 begin{eqnarray} p(theta|x)=frac{p(x|theta) p(theta)}{p(x)} end{eqnarray} ちなみに、(p(x))は以下のようにしておくとわかりやすい。 同時確率と周辺確率の関係。表を書いて縦、横がクロスするところが同時確率だとして、 縦、横いずれかの方向に同時確率を足し合わせる操作にあたるらしい。 なにか確率変数が独立でなければならない、というのは気にしない。 begin{eqnarray} p(x) = int p(x,theta) dtheta end{eqnarray} なので以下みたいに書き直せる。最後の比例のところは...。 左辺は事後確率分布、右辺は尤度と事前確率分布の積!!。 begin{eqnarray} p(theta|x) &=& frac{p(x|theta) p(theta)}{p(x)} \\ &=& frac{p(x|theta) p(theta)}{int p(x,theta) dtheta} \\ &propto& p(x|theta) p(theta) end{eqnarray} (p(theta))は確率変数(theta)の確率分布。尤度(p(x|theta))は(theta)に対して定数。(x)に対して変数。 ということで、右辺は確率分布(p(theta))を尤度(p(x|theta))を使って変形した確率分布。 で、左辺の(p(theta|x))は右辺の(p(x|theta) p(theta))を定数倍した確率分布。 データを観測していない状態で立てた(p(theta))があって、 観測したデータを使って求めた尤度(p(x|theta))が得られたことで、 左辺の(p(theta|x))が得られた、という状況。 (p(theta|x))は確率変数(theta)の確率分布なので、 最尤推定とベイズの定理を俯瞰してみると、最尤推定が点推定である一方で、 ベイズの定理では確率分布が得られるという具合で異なる。 (観測値が極端なデータだったとき、最尤推定は極端な推定結果が得られるだけだけれども、 ベイズの定理で得られる事後確率分布は確率分布なので様子がわかる..??) 事後確率分布を最大化する(theta_{MAP})を求めるのがMAP推定。(点推定) begin{eqnarray} hat{theta_{MAP}} = argmax_{theta} p(theta|x) end{eqnarray} 尤度をかけて得られた事後分布と同じ形になる便利な分布があって、 観測データ達の分布と対応して決まっている(共役事前分布)。 ベルヌイ分布の共役事前分布はベータ分布。

以前、\"正規分布に従う確率変数の二乗和はカイ二乗分布に従うことの証明\"という記事を書いた。 記事タイトルの通り、正規分布に従う確率変数の二乗和はカイ二乗分布に従う。 [clink url=\"https://ikuty.com/2018/08/01/chi-square-distribution\"] 実際にデータを生成して確かめてみる。 まずは、scipi.stats.chi2.pdfを使って各自由度と対応する確率密度関数を書いてみる。 \"pdf\" ってPortableDocumentFormatではなく、ProbabilityDensityFunction(確率密度関数)。 import numpy as np import matplotlib.pyplot as plt from scipy import stats # 0から8まで1000個のデータを等間隔で生成する x = np.linspace(0, 8, 1000) fig, ax = plt.subplots(1,1) linestyles = [\':\', \'--\', \'-.\', \'-\'] deg_of_freedom = [1, 2, 3, 4] for k in deg_of_freedom: linestyle = linestyles[k-1] ax.plot(x, stats.chi2.pdf(x, k), linestyle=linestyle, label=r\'$k=%i$\' % k) plt.xlim(0, 8) plt.ylim(0, 1.0) plt.legend() plt.show() 定義に従って各自由度ごとに2乗和を足し合わせてヒストグラムを作ってみる。 (ループのリストが自由度) cum = 0 for i in [1,2,3,4]: # 標準正規分布に従う乱数を1000個生成 x = np.random.normal(0, 1, 1000) # 2乗 x2 = x**2 cum += x2 plt.figure(figsize=(7,5)) plt.title(\"chi2 distribution.[k=1]\") plt.hist(cum, 80) 全部k=1ってなってしまった...。左上、右上、左下、右下の順に1,2,3,4。 頻度がこうなので、各階級の頻度の相対値を考えると上記の確率密度関数の形になりそう。 k=1,2が大きくことなるのがわかるし、3,4と進むにつれて変化が少なくなる。

[mathjax] 球面集中現象理解のための数学シリーズ第2弾。 前の記事でデカルト座標->極座標の変換から体積要素の積分により3次元球体の体積を導出してみました。 極座標の3変数(r,phi_1,phi_2)について定積分を計算していくと(frac{4pi r^3}{3})が出てきます! [clink url=\"https://ikuty.com/2019/08/02/concentration_on_the_sphere_1/\"] より高次元の超球体の体積がわかると、超高次元のときに体積が超球体の表面に集中する様が 論理の飛躍なく理解できるらしいので、今回は高次元の超球体の体積を導出するやつを流してみます。 (n)次元空間において以下を満たす点の集合が半径(r)の(n)次元球体です。 (x,y,z...)とやっていくとアルファベットが尽きるので(x_1,x_2,cdots,x_n)とします。 (3次元だとすると、半径が(sqrt{x_1^2+x_2^2+x_3^2}le r)の全ての点。) begin{eqnarray} x_1^2+x_2^2+cdots+x_n^2 le r^2 end{eqnarray} この集合を全部足したものが体積になります。 例えば3次元のとき、前回の記事の通り以下のようになります。ここではデカルト座標(x,y,z)。 極座標変換することで計算ができます。 begin{eqnarray} int_{x^2+y^2+z^2 lt r^2} U(x,y,z) dV &=& int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{r} U(r,phi_1,phi_2)r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 \\ &=& frac{4}{3} pi r^3 end{eqnarray} ちなみに2次元球体の体積(つまり円の面積)は当たり前ですが(pi r^2)です。 1次元球体の体積(つまり直線の長さ)は(r)です。 なんとなく(n)次元だと(r)の次数が(n)になりそうですが、 実際そうで(時間がないので公式として使ってしまいます)、 (n)次元の超球体の体積は(r^n)に比例します。 (Vr^n)のようにしておきます。 (n)次元の超球体を輪切りにすることを考えます。 球をどうやって輪切りにしても断面図は真円ですので...。 3次元の球体を輪切りにすると2次元の平面が現れる、というイメージです。 (n)次元の超球面を輪切りにし(n-1)次元の平面を作成し、 (n-1)次の平面と微小距離(Delta r)をかけて(n)次元の直方体を作ります。 微小距離(Delta r)を限りなくゼロに近づけることで体積要素とします。 体積要素は詰まるところ輪切りですが、輪切りなのに(n)次元なのです。 体積要素を全範囲で積分することで体積を求めます。 さて、(x_1^2 + x_2^2 + cdots + x_n^2 = 1)が半径1の単位球です。 これを輪切りにします。(n)次の変数を定数にします。(x_n=t)。 つまり、輪切りは(x_1^2 + x_2^2 + cdot + x_{n-1}^2 = 1-t)。 輪切りの半径は(sqrt{1-t})。(t=0)であれば丁度球体を真っ二つにする感じです。 (t=1)であれば球体のキワの限界で面積がゼロ。 (t)で輪切りにしたあと、(t)を微小距離(t+Delta t)だけ増やします。 その体積は(V_{n-1}(sqrt{(1-t)^{n-1}}) Delta t)。 (Delta t)を限りなくゼロに近づけると真の輪切りに近づき、 それを(r)の全範囲で定積分します。 begin{eqnarray} V_n=int_{-1}^{1} V_{n-1}sqrt{(1-t)^{n-1}} dt end{eqnarray} 球体は真ん中で対照なので、 begin{eqnarray} V_n =2 int_{0}^{1} V_{n-1}sqrt{(1-t)^{n-1}} dt end{eqnarray} 計算してから積分するのと、積分してから計算するのが同じなので、 begin{eqnarray} V_n=2 V_{n-1} int_{0}^{1} sqrt{(1-t)^{n-1}} dt end{eqnarray} 超絶懐かしい置換積分を使うと(int f(x) dx = int f(g(t)) frac{dx}{dt} dt)が出来る。 (t=sintheta)とすると、 begin{eqnarray} V_n= 2 V_{n-1}(int_{0}^{frac{pi}{2}} cos^{n-1}theta cos theta d theta ) end{eqnarray} これどうやって積分するんだよ...、と呆然としてしまうけれども、 公式があるようです。時間がないので公式を使います! 階乗が2つ並んでいるのは1つ飛ばしで階乗をする2重階乗というらしい。 begin{eqnarray} int_{0}^{frac{pi}{2}} sin^n x dx &=& int_{0}^{frac{pi}{2}} sin^{n-1}x sin x dx \\ &=& frac{pi}{2} frac{(n-1)!!}{n!!} end{eqnarray} ちなみに(cos^n x)の積分は(sin^n x)の積分から導かれて以下のようになる。 begin{eqnarray} int_{0}^{frac{pi}{2}} sin^n x dx = int_{0}^{frac{pi}{2}} cos^n x dx = begin{cases} frac{(n-1)!!}{n!!} & nが奇数 \\ frac{pi}{2} frac{(n-1)!!}{n!!} & nが偶数 end{cases} end{eqnarray} 体積の話に戻ると、 begin{eqnarray} V_n = begin{cases} 2V_{n-1}frac{(n-1) (n-3) (n-5) cdots 2}{n (n-2) cdots 3} & nが奇数 \\ pi V_{n-1} frac{(n-1)(n-3)(n-5)cdots 3}{n (n-2) (n-4) cdots 2} & nが偶数 end{cases} end{eqnarray} 偶数のときと奇数のときで別になってしまっているのを1つにしたい。 偶数のとき(n=2k)、奇数のとき(n=2k-1)とおいて、それが等しいとする。 begin{eqnarray} V_{2k-1} &=& 2V_{2k-2} frac{(2k-2)(2k-4)(2k-6)cdots 2}{(2k-1) (2k-3) cdots 3 } \\ V_{2k} &=& pi V_{2k-1} frac{(2k-1)(2k-3)(2k-5)cdots 2}{2k (2k-2)(2k-4) cdots 2} end{eqnarray} 2つの式で(V_{2k-1})が現れるので、それで等式を立てる。 begin{eqnarray} 2V_{2k-2} frac{(2k-2)(2k-4)(2k-6)cdots 2}{(2k-1) (2k-3) cdots 3 } = frac{V_{2k}}{pi} frac{2k (2k-2)(2k-4) cdots 2}{(2k-1)(2k-3)(2k-5)cdots 2} end{eqnarray} ガンガン約分する。 begin{eqnarray} V_{2k} = 2pi V_{2k-2} frac{1}{2k} = frac{pi V_{2k-2}}{k} end{eqnarray} (V_{2k})と(V_{2k-2})の漸化式になっていて、(V_{2k-2})が(V_2)になるまで再帰的に計算する。 begin{eqnarray} V_{2k} = frac{pi^{k-1}}{k!} V_2 = frac{pi^k }{k!} end{eqnarray} 階乗の一般化(ガンマ関数,(n! = Gamma{(n+1)}))を使って書くと、 begin{eqnarray} V_n = frac{pi^{frac{n}{2}}}{Gamma{(frac{pi}{2}+1)}} end{eqnarray} (n)次の超球面の体積が出ました...。 次回、これを使って超高次元でメロンパンの皮が厚いのを示します..。

[mathjax] 球面集中現象を理解するために記憶にないことが多すぎるので 理解に必要そうなことを少しずつ復習していきます。 まず極座標の直行座標変換。極座標を使って球の体積を求めてみます。 次の記事で球の体積を求める時に使う微小体積がヤコビアンになることを確認します。 極座標の直行座標変換 2次元のとき、極座標((r),(phi)) を直行座標で表すと、 begin{eqnarray} x &=& r cos(phi) \\ y &=& r sin(phi) end{eqnarray} 3次元のとき、極座標((r),(phi_1),(phi_2))を直行座標で表すと、 begin{eqnarray} x &=& r cos(phi_1) \\ y &=& r sin(phi_1) cos(phi_2) \\ z &=& r sin(phi_1) sin(phi_2) end{eqnarray} ここまではイメージできるのだけれども、(N)に飛ばしたときにどうなるのか。 N次元のとき、極座標((r),(phi_1),(phi_2),(cdots),(phi_{N-1}))を直行座標で表すと 以下のようになるらしい。 begin{eqnarray} x_1 &=& r cos(phi_1) \\ x_2 &=& r sin(phi_1) cos(phi_2) \\ x_3 &=& r sin(phi_1) sin(phi_2) cos(phi_3) \\ vdots \\ x_{N-1} &=& r sin (phi_1) cdots sin(phi_{N-2}) cos(phi_{N-1}) \\ x_{N} &=& r sin (phi_1) cdots sin(phi_{N-2}) sin(phi_{N-1}) end{eqnarray} この後使わないけれど、逆変換は以下の通りになるらしい。 begin{eqnarray} r &=& sqrt{x_1^2 + cdots + x_{N-1}^2 + x_N^2} \\ phi_1 &=& arccosfrac{x_1}{sqrt{x_1^2 + cdots + x_{N-1}^2 + x_N^2}} \\ phi_2 &=& arccosfrac{x_2}{sqrt{x_2^2 + cdots + x_{N-1}^2 + x_N^2}} \\ vdots \\ phi_{N-2} &=& arccos frac{x_{N-2}}{sqrt{x_{N-2}^2 + x_{N-1}^2 + x_N^2}} \\ phi_{N-1} &=& begin{cases} arccos frac{x_{N-1}}{sqrt{x_{N-1}^2+X_N^2}} & X_N ge 0 \\ -arccos frac{x_{N-1}}{sqrt{x_{N-1}^2+X_N^2}} & X_N lt 0 end{cases} end{eqnarray} 3次元極座標系における体積計算の例 極座標((r),(phi_1),(phi_2))において(r)を微小量(Delta r)だけ増やす。 (2pi r)は半径を(r)とする円の直径だということを考慮すると、 (phi_1)から微小角(Delta phi_1)増やしたときの(r Delta phi_1)は円弧の一部。 (Delta phi_1)を限りなく小さくすると(r Delta phi_1)は直線と考えられる。 (2pi rsin(phi_1) )が半径を(rsin(phi_1))とする円の直径と考えると、 (phi_2)から微小角(Delta phi_2)増やしたときの(r sin(phi_1) Delta phi_2)は円弧の一部。 (Delta phi_2)を限りなく小さくすると(r sin(phi_1) Delta phi_2)は直線と考えられる。 (Delta r, Delta phi_1, Delta phi_2)が限りなくゼロに近くときに 上記のそれぞれの直線を辺とする直方体ができる。その体積を(dV)とすると、 begin{eqnarray} dV &=& Delta r cdot r Delta phi_1 cdot r sin(phi_1) Delta phi_2 \\ &=& r^2 sin(phi_1) Delta r Delta phi_1 Delta phi_2 end{eqnarray} 直交座標(x,y,z)の位置を(U(x,y,z))とする。 直交座標系における球の体積はざっくり以下で表せる。 begin{eqnarray} int U(x,y,z) dV end{eqnarray} 直行座標の極座標変換は前述の通り以下。この座標を(U(r,phi_1,phi_2))とする。 begin{eqnarray} x &=& r cos(phi_1) \\ y &=& r sin(phi_1) cos(phi_2) \\ z &=& r sin(phi_1) sin(phi_2) end{eqnarray} 極座標系において以下だから、 begin{eqnarray} 0 le r le 1 \\ 0 le phi_1 le pi\\ 0 le phi_2 2pi end{eqnarray} それぞれについて積分すると体積になる。 begin{eqnarray} int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{1} U(r,phi_1,phi_2)r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 end{eqnarray} つまり以下。 begin{eqnarray} int U(x,y,z) dV = int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{1} U(r,phi_1,phi_2)r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 end{eqnarray} ここで( U(r,phi_1,phi_2) =1 )として積分範囲を(r)とすると球の体積の公式になる。 begin{eqnarray} int U(x,y,z) dV &=& int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{r} r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 \\ &=& int_{0}^{2pi} int_{0}^{pi} Bigl( int_{0}^{r} frac{sin(phi_1)}{3} bigl[ r^3 bigr]_{0}^{r} Bigr) dphi_1 dphi_2 \\ &=& int_{0}^{2pi} int_{0}^{pi} frac{r^3 sin(phi_1)}{3} dphi_1 dphi_2 \\ &=& int_{0}^{2pi} -frac{r^3}{3} bigl[ cos(phi_1) bigr]_{0}^{pi} dphi_2 \\ &=& int_{0}^{2pi} frac{2r^3}{3} dphi_2 \\ &=& frac{2r^3}{3} bigl[ phi_2bigr]_{0}^{2pi} \\ &=& frac{4pi r^3}{3} end{eqnarray} 次回、微小体積とヤコビアンが等しい理由を書いてみます。

[mathjax] sihoutteって何て読むのか...と思うけども\"シルエット\"だそう。フランス語源。 ポートレート写真を背景白、顔を黒に減色した例の\"シルエット\"。\"輪郭\"みたいな。 データ達を複数のクラスタに分割したとして、各々のデータがそのクラスタに存在することの 収まりの良さを表すことができる。 本来(相対的に)他のクラスタに属しているべきデータ達の割合が分かったり、 クラスタ境界で(相対的に)どちらのクラスタに属しても良さそうなデータ達の割合が分かったりする。 sklearnは教師なしクラスタリングまでで、それをグラフ化しないといけないのだけど、 matplotlibに対応するのがなく、\"線を書く\"みたいなコマンドを並べて作っていかないといけない様子。 yellowbrickというパッケージを使うとそれもやってくれる。 Yellowbrick is a suite of visual diagnostic tools called “Visualizers” that extend the Scikit-Learn API to allow human steering of the model selection process. In a nutshell, Yellowbrick combines scikit-learn with matplotlib in the best tradition of the scikit-learn documentation, but to produce visualizations for your models! sklearnとmatplotloibだけで出力するバージョンと、yellowbrickで出力するバージョンの 両方を試してみた(前者はクラスタのラベルを出力できずsihoutteグラフとの対応関係が理解できずに 中途半端で終了)。 sihoutte係数 全データ達は(x_i)。データ(x_i)がクラスタ(A)に属し、最近傍にクラスタ(B)があるという状況。 (A)のクラスタ中心は(x_A)、(B)のクラスタ中心は(x_B)。 \"最近傍のクラスタ中心との距離の平均\"から\"自分のクラスタ中心との距離の平均\"を引いた値。 わざわざ1行に全ての変数が出てくるように書いてみる。 begin{eqnarray} s_i = frac{sum_{i=1}^{n}|x_i-x_B|/n - sum_{i=1}^{n}|x_i-x_A|/n }{max bigl{sum_{i=1}^{n}|x_i-x_A|/n,sum_{i=1}^{n}|x_i-x_B|/n bigr}} end{eqnarray} データが(n)個あるので、sihoutte係数も(n)個できる。 自分が属しているクラスタ(A)よりも、隣のクラスタ(B)の方が居心地が良いと 分子はマイナスになる。 今属しているクラスタでも隣のクラスタでも、どちらもさほど居心地が変わらないとゼロ近辺になる。 大きければ迷いなく今のクラスタで良いことを表せる。 sklearnとmatplotlibだけでsihoutte係数 乱数で作った偽データに6-meansをかけた図。 どうやってもクラスタ中心のラベルが取り出せない!!..(スルー..) import numpy as np import matplotlib.pyplot as plt import seaborn as sns mu = [[0,0], [20,20], [50,50], [40,30], [40,10], [20,40]] sigma = [ [[30,20],[20,50]], [[20,30],[10,20]], [[60,40],[20,20]], [[60,20],[20,60]] ,[[30,10],[10,30]],[[50,20],[20,50]] ] points = 100 clusters = [] for index in range(len(mu)): cluster = np.random.multivariate_normal(mu[index], sigma[index], points) dig = np.full((points,1),index+1, dtype=int) cluster = np.hstack((cluster,dig)) clusters = np.r_[clusters,cluster] if len(clusters) > 0 else cluster plt.scatter(x=clusters[:,0], y=clusters[:,1],c=clusters[:,2]) from sklearn.cluster import KMeans model = KMeans(n_clusters=6, init=\'random\',max_iter=10) y_km = model.fit_predict(clusters[:,:2]) labels = model.labels_ centers = model.cluster_centers_ centers fig = plt.figure() ax1 = fig.add_subplot(1,1,1) ax1.scatter(x=clusters[:,0], y=clusters[:,1],c=labels) ax2 = fig.add_subplot(1,1,1) ax2.scatter(x=centers[:,0], y=centers[:,1], alpha=0.5,s=600,c=\"pink\",linewidth=2,edgecolors=\"red\") sklearnとmatplotlibだけでsihoutte係数のグラフを出してみる。 \"線を書く\"みたいなコマンドを並べて規格通りの図を作るのか...。 from sklearn.metrics import silhouette_samples # cluster数 num_clusters=6 #全データのsilhouette係数を取得 silhouette_vals = silhouette_samples(clusters[:,:2],y_km,metric=\'euclidean\') cluster_labels = np.unique(y_km) min,max = 0,0 yticks = [] for i,c in enumerate(cluster_labels): c_silhouette_vals = silhouette_vals[y_km == c] c_silhouette_vals.sort() max += len(c_silhouette_vals) plt.barh( range(min,max), c_silhouette_vals, height=1.0, edgecolor=\'none\' ) yticks.append((min+max)/2.) min += len(c_silhouette_vals) avg = np.mean(silhouette_vals) plt.axvline(avg, color=\'red\', linestyle=\"--\") plt.yticks(yticks, cluster_labels + 1) plt.ylabel(\'Cluster\') plt.xlabel(\'Silhouette coefficient\') plt.show() 全体的に茶色のクラスタのsilhouette係数が低め。 残念ながら、どのクラスタと対応するのか出力できず何の考察も出来ず...。 yellowbrickでsilhouette係数を出力 無茶苦茶簡単に出せる。 from yellowbrick.cluster import SilhouetteVisualizer sv = SilhouetteVisualizer(model) sv.fit(clusters[:,:2]) sv.poof() 出てきた図。自作したものと全然合っていないように見える.. 自作したものとラベルが合っていないだけだと信じたい。 似た形の塊があるので.. 以上、完全な失敗だけれども一旦終了...