２つの株を平均的に購入するための株数を計算する
問題
A君は月毎に一定額を貯蓄をしようと考えましたが、銀行の金利が低すぎるのでインデックスファンド（ETF)への積み立てを考えました。ファンドAとファンドBにいままでの額と合わせて評価額が同額になるように毎月割り当てようとしましたが、10口単位での購入のため決まった額範囲内でバランスをとる必要があります。毎月の時価に合わせてそれぞれのファンド何口ずつかう必要があるでしょうか。 <<
という悩みをプログラムで解決してみます
積み立てをおこなっているのですが、国内インデックスと海外インデックスのバランスを取りながら毎月購入をしたいんですよね。購入時期になった時の悩みが10口単位の口座を何ずつ買えばよいのか、ということになります。解決するためのプログラムを書いてみました。Git Hubにコード上げてます。
解説
まず本日の価格を設定する必要があります。関数で適当に設定しますが、一応ランダムに変更するように設定しました。ファンドは2つあるので数値で識別します。
def stock_value(stk): if stk == 0: value = 1400 elif stk == 1: value = 1000 err = random.randrange(-100, 100) return value + err
あとは初期設定。現在持っているそれぞれのストック数、価格などを適当に設定します。
# Current stocks you have purchased_stock = [90, 110] # range to purchase in next time purchase_range = [50000, 60000] # amount of unit to buy in next time buy_balance = [0, 0] # List for buy units balances = [] # Dispersion vers = [] # current value for each account stk_value = [stock_value(0), stock_value(1)]
10口単位でそれぞれの株数を上げながら、購入する価格のレンジ内に収まる株数をチェックします。総当たりです。
total_value = sum([buy_balance[i]*stk_value[i] for i in range(2)]) while total_value < purchase_range[1]: while total_value < purchase_range[1]: # If the balance is in purchase range, # culculate dispersion and add lists if purchase_range[0] < total_value : balances.append(list(buy_balance)) values = sum([(purchased_stock[i]+buy_balance[i]) * stk_value[i] for i in range(2)]) ver = sum([((purchased_stock[i]+buy_balance[i])*stk_value[i] - values/2)**2 for i in range(2)]) print(buy_balance, stk_value, total_value, values) vers.append(ver) buy_balance[0] += 10 total_value = sum([buy_balance[i]*stk_value[i] for i in range(2)]) buy_balance[0] = 0 buy_balance[1] += 10 total_value = sum([buy_balance[i]*stk_value[i] for i in range(2)])
それぞれの組み合わせの中で分散値が一番低いものを選びます。それが、総額のバランスがとれた今回の組み合わせの口数を示す番号となります。
# Checking minimum dispersion in the list idx = vers.index(min(vers))
結果
BALANCED PLANが今回買うべき株数になります。
======== RESULT ======== CURRENT STOCK VAL: [1440, 943] CURRENT STOCK : [90, 110] CURRENT PURCHASED: [129600, 103730] CURRENT VALUE : [1440, 943] PURCHASE RANGE : [50000, 60000] BALANCED PLAN : [10, 40] AFTER PURCHAE : [144000, 141450]
漫画「一丁目のトラ吉」
ブックマークコメントでも書いたのですが、「一丁目のとら吉」という漫画が実家に置いてあってよく読んでました。誰がどのようにして買ったのか全く不明だったのですが。 この増田をみてついつい思い出して検索したところ、電子書籍でもあるんですね。本当にいい時代。 5巻完結、そういえば途中までしかなくて続きを読んでいなかったこととも思い出してついつい一括購入しちゃいました。
![[まとめ買い] 一丁目のトラ吉](https://images-fe.ssl-images-amazon.com/images/I/B17aubpzaAS._SL160_.png "[まとめ買い] 一丁目のトラ吉")
- 作者: セツコ・山田
- メディア: Kindle版
- この商品を含むブログを見る
主人公は猫のこの作品。さらに特徴として猫同士は言葉によるコミュニケーションがとれるのだが、人間の言葉は象形文字のような言語で、「なにか言ってるんだろうけど、内容は詳しく理解できない」様子が表現されています。ネコ語、人間語なる解釈もでてきます。
この物語は飼い主家族や妹猫のチビ、犬のモジャとの「家族・家庭」、さらには飼い猫、野良をふくめたご近所さんの猫たちとの「外」とのコミュニティを舞台として、人間社会に適応して生活を送る猫たち物語をトラ吉の目線からすすめられていきます。何気ない人間の日常が猫目線になると途端に滑稽に見えてくるような描写や、捨て猫、責任感が薄い飼い主など猫を取り巻く日常の中の社会問題が当事者？の視点から書かれているので実はよりリアリティが浮きだってくるように見えます。
私が特に好きな話はシロクロという猫の回です。 ごはんだけ特定の家からいつももらって普段は野良のように生活をしている半野良のシロクロは、その飼い主家族が自分をおいて突然引っ越しをしてしまったことを知ります。
それでも、もぬけの殻となった家の前を離れようとしないシロクロ。まわりの猫も彼を諭すのですが頑に居座り続けます。そんななかその家に別の家族が引っ越してきます。シロクロとその新しい家族とのやりとりがトラ吉目線で進められます。
作者のセツコ・山田さんは自分の家族や飼い猫をモデルにして実際に起きたことをベースに猫目線で描いているのかと思われます。猫を愛し、よく観察していなくてはこんなリアルな話はなかなか書けないのではないでしょうか。
猫好きにこそぜひ読んでいただきたい作品となります。
反復性肩脱臼の手術から1年経過
早いもので反復性肩脱臼からおよそ一年たちました。
このブログでは機械学習系の記事をメインで書いていますが、いまだに半分以上のアクセスがこちらの記事のようです。同じ悩みを抱える人は意外と多く、どのような手術なのか気になって人も多いみたいですね。こういうレポートって大事なんだなとアクセスログを見るたびに思います。そこで一年後どうなっているかの記録と、思い出せる限りの1年の経過も残しておきたいと思います。
現在、1年後の状況
当然ですが、不自由なく普通に生活ができています。手術した左肩はきちんと後ろまでまわせますし、脱臼することはないです。可動域は右と比べて少ないのは否めないと思います。後ろに回すと引っ張られる感じはありますが、人口のひもで固定しているので仕方ないのかと思います。あとリハビリを頑張ればもっと行くのかもしれませんが、私生活でその必要は感じられないので私個人としては十分といえます。
一年の経過
- 術後すぐ：術部の回復をおこなうため、サポータによる固定で過ごします。利き手ではなかったので、まだよかったですが、それでも不自由でした。寝るときもつけたままでねるので、寝返りもできないし、ずっと吊っている左肩の肩こりはとてもひどかったのを覚えています。リハビリも週に2回です。
- 術後一か月経過：サポーターは外せますが、可動域はとても狭いです。前から上に手を上げようとしても最初は90度まであがりませんでした。生活にも影響がでます。左手を使って胴体から離して行わなくてはいけない作業は大変です。高いものを取るとかは右手でもできますが、何気に靴ひもを結ぶことは両手を使わなくてはならず、さらに胴体から離れてやる作業なので難しいです。リハビリは引き続き週二回です。
- 術後三か月経過：だんだんと可動域が広がっていきますがまだ制限がでます。前から上に腕を上げる運動は180度までもう少しのところで止まってしまい、なかなかそこから改善せずもどかしい思いをしました。診断の時もリハビリが芳しくないということでもし次回の診断までによくなければ筋肉を柔らかくする注射をうつといわれました。結果として一か月後の診断では合格点がでてその注射はうたなくてもよくなりましたが。それを受けてリハビリの回数も週1回になりました。
- 術後六か月経過：運動が解禁されて、ジム通いもスタートしました。可動域も実生活に問題がない程度まで動くようになりました。そこから数か月してリハビリも終わりとなりました。
というわけでリハビリが不安になった時もありましたが、結果として日常生活にはなにも問題ないレベルまで回復しました。脱臼するかもしれないと可動域をセーブしていた術前が嘘のように特に何も気にせずに生活をしています。半年以上のリハビリはありましたが、今後の人生や老後のことも考えるとやってよかったとはっきりと言えます。お悩みの方は半年以上はすこし不自由ですが、なるべく早いうちに決断することをおすすめします。
<SVM> サポートベクタマシン自主勉強（１）：方針
ちょっと間が空いてしまいましたが機械学習のエントリとなります。前回最後にちらっとお話ししたのですが、サポートベクタマシン（SVM）に足を突っ込んでみたいと思っています。
SVMのコーディングが自力でできない問題
ただ、このSVMは何回か挑戦して失敗している経緯があります。前回のようにはじパタ本だけでは私のレベルでクリアができない代物となります。そこではじパタと合わせて、「異常検知と変化検知」の内容と相互補完してみています。
- 作者: 井手剛,杉山将
- 出版社/メーカー: 講談社
- 発売日: 2015/08/08
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (2件) を見る
そこでやっと、SVMの概要について理解した。。。。つもりでした。しかしながらいざこれをスクリプトに落とそうとするとどう書いていいのかわからない。つまりロジックがわかったつもりでも、実行に移すとなるとどうしていいのかわからなくなってしまいました。
頼るべきはGoogle先生と集合知
しかしながらもうすでに歴史のあるこの分野では検索してしまえばいくらでもコードが出てきます。今回は以下のサイトさまのコードを拝借してやはりアヤメデータを使った分類を行ってみました。
こちらのコードは の線形で分離をして
を1, -1とわけています。アヤメデータもとりあえず、3種のうちの最初の1種を1, その他を-1としてインプットするコードに改造をしてみました。
その結果が以下となります。
おぉ、なんとなく分類する線が学習を経るごとに赤と緑の間にはいっていきますね。赤が1としたアヤメデータなのでうまくいっているようです。
ここからスタート
SVMについてまとめることはたくさんありそうです。ラグランジュ乗数法、正則化、双対問題、KKT条件、カーネルトリック。数学的要素が多くどこまでカバーできるかわかりませんが、まず次回はこのスクリプトについて解いていくことから初めてみようかと思います。正直きちんと着地できるかまだ不明瞭ではあります。
コード
Jupyter Notebookで実行しています。
import sys import numpy as np from sklearn import datasets import itertools import matplotlib.pyplot as plt #inline pyplot from matplotlib import animation from IPython.display import HTML iris = datasets.load_iris() features = iris.data feature_names = iris.feature_names targets = iris.target class SVM(): def __init__(self, v1=0, v2=2, t=0): self.X = np.array([features.T[v1], features.T[v2]]).T self.T = np.array([1 if i==t else -1 for i in targets]) self.N = self.X.shape[0] self.d = self.X.shape[1] self.alpha = np.zeros(self.N) self.beta = 1.0 self.eta_al = 0.0001 self.eta_be = 0.1 self.v = [v1, v2] def test(self): print(self.X.shape) print(self.T) print(self.N) print(self.alpha) def plot_map(self, w, b): plt.subplot(1, 1, 1) seq = np.arange(4.0, 8.0, 0.02) plt.xlabel(feature_names[self.v[0]]) plt.ylabel(feature_names[self.v[1]]) plt.xlim(4, 8) plt.ylim(0, 8) plt.plot(seq, -(w[0] * seq + b) / w[1], 'k-') #plt.autoscale() plt.grid() for t, marker, c in zip(range(3), '>ox', 'rgb'): #print(self.X[self.T==1], 0) #print(self.T) plt.scatter( features[targets == t, self.v[0]], features[targets == t, self.v[1]], marker=marker, c=c, ) #print(linear_pro(w, np.array([1.,1.,1.]))) plt.show() def rtin(self, show=False): k_w = np.array([-100.0, -100.0]) k_b = -100.0 value = [] for _itr in range(1000): for i in range(self.N): self.delta = 1 - (self.T[i] * self.X[i]).dot(self.alpha * self.T * self.X.T).sum() - self.beta * self.T[i] * self.alpha.dot(self.T) self.alpha[i] += self.eta_al * self.delta for i in range(self.N): self.beta += self.eta_be * self.alpha.dot(self.T) ** 2 / 2 index = self.alpha > 0 w = (self.alpha * self.T).T.dot(self.X) b = (self.T[index] - self.X[index].dot(w)).mean() value.append([w, b]) if (np.power(np.array([k_w - w]), 2) < 0.00001).all() and (k_b - b) ** 2 < 0.00001: print(_itr) break if _itr % 10 == 0 and show: self.plot_map(w, b) k_w = w k_b = b return value svm = SVM(t=0) _ = svm.rtin(True)
＜はじパタ＞ おまけ：混合ガウスモデルの検証
前回の記事より。 bython-chogo.hatenablog.com
EMアルゴリズムの実装を行ってみたのですが、アヤメの2次元データではこの混合ガウスモデル、EMアルゴリズムはうまく適用できないね、という結論になりました。2つの山が重なったピーナッツを横にクラスタリングしてほしいのに、縦に割ってしまうようで、初期値をランダムにするとほとんどうまくいかにようでした。
でも、本当でしょうか。私のコードになにか問題があって、うまくクラスタリングができなかった可能性を排除しきれません。いや、むしろこんな超適当コードを信じられないのです。そうだ、私のコードが間違っていてちゃんとEMアルゴリズムの実装ができればうまくいくのではないか。
そんな疑念についてはっきりさせるために積読してあった以下の書籍に書かれているコードを改造して改めて同じデータのプロットをしてみました。
Pythonで学ぶあたらしい統計学の教科書 (AI & TECHNOLOGY)
- 作者: 馬場真哉
- 出版社/メーカー: 翔泳社
- 発売日: 2018/04/19
- メディア: 単行本（ソフトカバー）
- この商品を含むブログを見る
その結果が、こ・れ・だ
やっぱり縦に割ってるじゃねーか！
というわけでEMアルゴの実力はほかのデータセットで確認をしたほうがよさそうですね。次回ははじパタの11章ですが、もしかしたらサポートベクタマシンのほうに浮気するかもしれません。
＜はじパタ＞ 10.4 確率モデルによるクラスタリング
はい、ここまで10章の教師なし学習のクラスタリングをやってまいりましたが、10章も最後となりました。とりあえず一章をきちんと走り抜けられてよかったです。 で、ここは前回までのK-平均法や融合法と違い、各要素がそれぞれのクラスに入る確率を求める確率モデルによるクラスタリングをおこないます。キーワードとしては、混合正規分布モデルとEMアルゴリズムとなるようです。わかりやすい説明はこちらのリンクをご参照、とりあえず自分なりにまとめてみます。
混合正規分布モデル
各クラスタの分布がそれぞれ正規分布にしたがうと仮定をしてその平均と分散が一番Fitする最尤の値を求めていくものになります。また、ある要素に関してそれぞれクラスタの割合を表現する混合比なる重みづけのためのパラメータもでてきてこの値もいじる必要があります。
EMアルゴリズム
所属しているクラスタを示す隠れ変数があるものとして、その変数の確率モデルの最尤推定値を求めるのに優れているのが、EMアルゴリズムとなるようです。E（Expectation）ステップで隠れ変数に対しての期待値を計算、M（Maximization）ステップでその隠れ変数の最大値を計算するようにします。
はじパタに限らず、いくつかの本にこのEMのステップで何が行われているかの数式や説明があるのですが、いかに示すサイトの説明が図として直感的に何をしているのかわかりやすく書かれていると思いました。こういうの貴重。
www.slideshare.net
ピーナッツを縦に割っちゃう問題
それで、コーディングを走らせてみたのですが、どうしても正解率が上がらないんです。実際にプロットをしてみると以下のように右上のピーナッツを縦にわるようなクラスタリングをしてしまうのです。
これは初期値にだいぶよるようですね。平均値のmuの値を意図的に変えてみたところ以下のようにピーナッツを横にわるような「理想的な」クラスタに収束することもあるようです。
初期値をランダムに置くようにして確認してみましたが、この理想的な形状に収まるのは10回に1～2回程度のようです。
合わないモデル？
初期値問題に頭を抱えて機械学習のお師匠様のところに相談したところ、「そもそも合わないモデルだったのでは」とのことです。考察としてはアヤメデータ以外でもこのEMアルゴリズムを試してみてうまく当てはまるかを確認してみたほうが良いかもしれないです。やるかはわからないですが。
ソースコード
いつも以上にとっちらかったソースになってしまっていると思います。まあ備忘録なので。
import numpy as np from sklearn import datasets import itertools import matplotlib.pyplot as plt #inline pyplot iris = datasets.load_iris() features = iris.data feature_names = iris.feature_names targets = iris.target from collections import Counter class Iris_PModel: def __init__(self): self.p = [] self.x = [] self.y = [] self.t = [] self.pt = [] self.num = [50, 50, 50] iris = datasets.load_iris() features = iris.data targets = iris.target def install(self, a=2, b=3): self.x = features.T[a] self.y = features.T[b] self.t = targets self.pi = [1/3, 1/3, 1/3] self.means_dists() def plot_map(self, df=False): p = self.pt if df: p = self.t for t, marker, c in zip(range(3), '>ox', 'rgb'): plt.scatter( features[p == t, 2], features[p == t, 3], marker=marker, c=c, ) plt.scatter( self.mu[t][0], self.mu[t][1], marker='s', c=c ) plt.xlim(1.9, 4.5) plt.ylim(0.9, 7.1) plt.xlabel(feature_names[2]) plt.ylabel(feature_names[3]) plt.autoscale() plt.grid() plt.show() def randomize(self): self.pt = np.array([np.random.randint(3) for i in range(sum(self.num))]) def norm(self, mu, sigma): n_mu = np.array(mu) n_sig = np.array([[sigma[0], sigma[2]],[sigma[2], sigma[1]]]) n_det = np.linalg.det(n_sig) n_inv = np.linalg.inv(n_sig) nv = [] for i in range(150): n_c = np.array([self.x[i], self.y[i]]) - n_mu nv.append((1/(2*np.pi*np.sqrt(n_det))) * np.exp(- n_c.dot(n_inv).dot(n_c).T/2)) return(nv) def means_dists(self, init=True): if init: self.randomize() x_m = sum(self.x)/150 y_m = sum(self.y)/150 self.nu = [[np.random.normal(x_m, 1), np.random.normal(y_m, 1)], [np.random.normal(x_m, 1), np.random.normal(y_m, 1)], [np.random.normal(x_m, 1), np.random.normal(y_m, 1)]] print("mu", self.mu) self.v = [[1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [1.0, 1.0, 0.0]] def m_step(self, gamma): n_k = [sum([gamma[i][j] for i in range(150)]) for j in range(3)] self.mu = [[sum([self.x[i]*gamma[i][j] for i in range(150)])/n_k[j], sum([self.y[i]*gamma[i][j] for i in range(150)])/n_k[j]] for j in range(3)] self.v = [[sum([(self.x[i]-self.mu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j], sum([(self.y[i]-self.mu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j], sum([(self.x[i]-self.mu[j][0])*(self.y[i]-self.mu[j][1])*gamma[i][j] for i in range(150)])/n_k[j]] for j in range(3)] self.pi = [n_k[i]/150 for i in range(3)] def e_step(self, init=False): if init: self.install(2, 3) nv = [self.norm( self.mu[i], self.v[i]) for i in range(3)] snv = [nv[0][i]*self.pi[0] + nv[1][i]*self.pi[1] + nv[2][i]*self.pi[2] for i in range(150)] prop = [[nv[0][i]*self.pi[0]/snv[i], nv[1][i]*self.pi[1]/snv[i], nv[2][i]*self.pi[2]/snv[i]] for i in range(150)] return prop def process(self): gamma = self.e_step(True) self.m_step(gamma) pre_result = [] for i in range(31): result = [] for g in gamma: result.append(g.index(max(g))) self.pt = np.array(result) max_int = [0, 0, 0] if pre_result == result: print(self.pi) self.plot_map() break pre_result = list(result) gamma = self.e_step(False) self.m_step(gamma) krs = Iris_PModel() krs.process()
＜はじパタ＞ おまけ：共分散行列のコーディングに苦戦した話
今、はじパタのコーディングをせっせとやってはブログに上げているという作業を続けていますが、いかんせん付け焼刃の知識のために基本的と思われるところで引っかかってしまうことも、多く、今後の備忘録的にも失敗とリカバリの記録を残してみようと思います。
現象：NaNエラーが発生
10.4 「確率モデルによるクラスタリング」をやっていくと、EMモデルの繰り返しによって数値がNaNとなってしまう現象が発生。ひとつNaNがでてくるとその数値を使った値もNaNとなるので真賀田博士もびっくりのすべてがNaNになってしまうのでした。以下がその時のログ。
MU: [[4.2591681450807695, 1.3885253632923946], [3.3363395791482486, 1.0264189358663134], [3.6851362069583895, 1.1830943936011733]] V: [[2.8113007747631147, 2.9874538055910858, 3.042873375542265], [0.51266119363098883, 0.56519089148146107, 0.59171106814975671], [1.1252472995735809, 1.2722812742580845, 1.296157580183261]] N_DET: -0.860447181357 N_DET: -0.0603705511147 N_DET: -0.24839340451 NV: nan nan nan SNV: nan nan nan nan nan
原因究明
NaNとなってしまう原因はどこにあるのかを突き止めてみると、ガウス分布の計算内で分散行列の行列式を出すのですが、そこで負の値になってしまっているようでした。その値の平方根を取るものだから当然実数にはならず、NaNとなってしまったのです。なるほど。 そもそものこの共分散行列をかくにんしたところ、
{\displaystyle \Sigma } \Sigma は、半正定値行列
分散共分散行列 - Wikipedia
半正定値行列というのは、行列式が０か正の値になる行列のこと。ということはそもそも負の値が出てくること時点でおかしいということです。 よし、原因はつかめた。あとは訂正をするだけだ。。。ここからが長かった。
行列式のコーディングのチェック
問題となるのは以下の分散を求める部分のコード。意図的にNumpyの行列を使わず、理解のため分散、共分散をそれぞれ計算してみようということで以下のリスト内包表記のコードとなっています。
self.v = [[sum([(self.x[i]-self.nu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)], [sum([(self.y[i]-self.nu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)], [sum([(self.x[i]-self.nu[j][0])*(self.y[i]-self.nu[j][1])*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)]]
最初のsumと次のsumが2次元の要素それぞれの分散値を計算、最後のsumが共分散を計算しています。なんども確認して問題なさそうだったのですが、その結果が以下のように行列式をとると負の値がでてきてしまう。なぜだ。
#分散1、分散2、共分散、行列式（ad-bc）
2.81130077476 2.98745380559 3.04287337554 -0.860447181357
0.512661193631 0.565190891481 0.59171106815 -0.0603705511147
1.12524729957 1.27228127426 1.29615758018 -0.24839340451
入力値とかも確認してみましたがそもそも共分散行列の行列式が負になること自体が異常なので、やはりこのコードがおかしいんですよね。 で、最終的に別の本に書いてあるEMアルゴリズムをコーディングしてみて、その行列要素を出力してみることにしました。すると。。。
V: # Result of ML Book's code
2.69263977682 0.640353547689 1.26287084435 0.129398664238
2.95622008252 0.425917289351 1.08691249439 0.0777264738108
2.85750222648 0.516449549181 1.19418340666 0.0496817279099
V: # Result of my code
2.69263977682 2.95622008252 2.85750222648 -0.205283191105
0.640353547689 0.425917289351 0.516449549181 0.00601751040866
1.26287084435 1.08691249439 1.19418340666 -0.0534439092193
おおぅ、行列の転置が起きてしまっているようですね。なるほど、リスト内包表記の組み合わせが正しくないようです。基本的な考え方は間違っていないかったですが、こうした順序の認識がまちがっていて、かつ訂正することも難しい（なにが間違っているか理解できない）状況に陥っていました。
訂正
これを踏まえて以下のように訂正をおこなったところ、出力からNaNが出てこなくなりました。
#Corrected self.v = [[sum([(self.x[i]-self.nu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j], sum([(self.y[i]-self.nu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j], sum([(self.x[i]-self.nu[j][0])*(self.y[i]-self.nu[j][1])*gamma[i][j] for i in range(150)])/n_k[j]] for j in range(3)] #Wrong self.v = [[sum([(self.x[i]-self.nu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)], [sum([(self.y[i]-self.nu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)], [sum([(self.x[i]-self.nu[j][0])*(self.y[i]-self.nu[j][1])*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)]]
反省
基礎知識がたりないですね。でもまぁ、もともと基礎知識が薄い中でこうしてコケることは織り込み済みで、コケながら覚えていくしかないかと思ってます。今回学んだことをまとめるといかになります。
- 共分散行列は半正定値行列。行列式が負の値になることはない。
- リスト内包表記は便利だけど順番を間違えやすい。転置状態に注意。
それでは次回はようやく、「確率モデルによるクラスタリング」にいけるかと思います。
