紙媒体で管理するとなくなりがちなのでブログで進捗などを管理することにしました
※殆どの記事は自分自身のためだけにかいています.他人に見せられるレベルには至っていません...

【Ubuntu】MeCabインストール

Ubuntu16.06LTSにMeCabとPython3用のMeCabのインストールを行った.

MeCabのインストール

1. Gitからcloneして展開する

$ git clone --depth 1 https://github.com/neologd/mecab-ipadic-neologd.git

2. mecab-ipadic-NEologdをインストールする

$ ./bin/install-mecab-ipadic-neologd -n -a

3. 実行する
自分の場合は/usr/local/lib/mecab/dic/mecab-ipadic-neologd/にあったため

$ macab -d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/

4. コマンドとして設定する
いちいちパスを宣言するのは面倒なので

$ alias mecab="mecab /usr/local/lib/mecab/dic/mecab-ipadic-neologd/"
Python3でMeCabを使う

1. 辞書のPATHを設定

sudo emacs /usr/local/etc/mecabrc
;dicdir = /usr/local/lib/mecab/dic/ipadic; これを消す
dicdir = /usr/local/lib/mecab/dic/mecab-ipadic-neologd;これを追加

2. GoogleDriveからPython3用のファイルを取得
https://drive.google.com/drive/folders/0B4y35FiV1wh7fjQ5SkJETEJEYzlqcUY4WUlpZmR4dDlJMWI5ZUlXN2xZN2s2b0pqT3hMbTQ
ここからMeCab-python3.*というファイルを取得し解凍

3. インストール
先ほど解凍したところに移って

$ python setup.py build
$ python setup.py install

4. 試しに使ってみる

import MeCab
mecab = MeCab.Tagger("-Ochasen")
print(mecab.parse("今日の天気は晴れです。"))

引用:http://blueskydb.blog.jp/archives/67055421.html


【参考】
http://blueskydb.blog.jp/archives/67055421.html
http://qiita.com/mojibakeo/items/2906162edf7045899974
http://qiita.com/knknkn1162/items/8c12f42dd167aae01c02
https://github.com/neologd/mecab-ipadic-neologd
https://pypi.python.org/pypi/mecab-python3
http://toriaezu-engineer.hatenablog.com/entry/2016/08/24/234145
http://qiita.com/yuichy/items/5c8178e5cc3711386b77

確率分布書き方

ポアソン分布の作り方
import scipy.stats
import numpy as np
import matplotlib.pyplot as plt
#figsize(12.5,4)

lambda_ = [1.5, 4.25]
colors = ["#348ABD","#A60628"]

a = np.arange(16)
plt.bar(a,scipy.stats.poisson.pmf(a,lambda_[0]),color=colors[0],label="$\lambda=%.1f$"%lambda_[0],alpha=0.60,edgecolor=colors[0],lw="3")
plt.bar(a,scipy.stats.poisson.pmf(a,lambda_[1]),color=colors[1],label="$\lambda=%.1f$"%lambda_[1],alpha=0.60,edgecolor=colors[1],lw="3")

plt.xticks(a+0.4,a)
plt.legend()
plt.ylabel("Probability of $k$")#kの確率
plt.xlabel("$k$")
plt.title("Probability mass function of a Poisson random variable, differing $\lambda$ values")
plt.show()

f:id:umashika5555:20170926155400p:plain

指数分布
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

a = np.linspace(0,4,100)
lambda_ = [0.5,1]
colors=["#348ABD","#A60628"]

for l,c in zip(lambda_,colors):
    plt.plot(a,scipy.stats.expon.pdf(a,scale=1./l),lw=3,color=c,label="$\lambda=%.1f$"%l)
    plt.fill_between(a,scipy.stats.expon.pdf(a,scale=1./l),color=c,alpha=.33)

plt.legend()
plt.ylabel("Probability density function at $z$")#zにおける確率密度関数
plt.xlabel("$z$")
plt.ylim(0,1.2)
plt.title("Probability density function of an exponential random variable, differing $\lambda$ values")
plt.savefig("exponential_distribution.png")

f:id:umashika5555:20170926161913p:plain

ベイズ推論イントロ

Pythonで体験するベイズ推論

Pythonで体験するベイズ推論 PyMCによるMCMC入門

Pythonで体験するベイズ推論 PyMCによるMCMC入門

  • 作者: キャメロンデビッドソン=ピロン,玉木徹
  • 出版社/メーカー: 森北出版
  • 発売日: 2017/04/06
  • メディア: 単行本(ソフトカバー)
  • この商品を含むブログを見る
から, ベイズ推論についての紹介

世界に司書と農家しか職業選択が無いときに司書と農家の人口比は1:20である.
Steveが内向的でおとなしい性格だという情報を得たときに, Steveは司書である確率はどのくらいなのであるのだろうか?

 P(A)をSteveが司書である確率とする.
このとき求めたいのはSteveが内向的な性格であるとわかった上で, Steveが司書であるという確率 P(A|X)である.
ベイズの定理より
 P(A|X) = \frac{P(X|A)P(A)}{P(X)}
 P(X|A)は司書であるときに内向的な性格であるという確率である. これは0.95とする.
 P(A)は司書である確率なので \frac{1}{21}である.
 P(X) = P(X かつ A) + P(X かつ not A) = P(X|A)P(A)+P(X|not A)P(not A)より,
 P(X|not A)は農家であるときに内向的な性格である確率で0.5とする.
すると P(X) = 0.95*\frac{1}{21} + 0.5*\frac{20}{21}より P(A|X) = 0.087となる.
これを図にすると下のようになる.
f:id:umashika5555:20170925165930p:plain
「Steveは内向的な性格だ」という事前情報を得たときの「Steveが司書である」という事後確率は少し高くなった.
コードは以下の通りである.

from IPython.core.pylabtools import figsize
import numpy as np
import matplotlib.pyplot as plt
figsize(12.5,4)
colors = ["#348ABD","#A60628"]
pr_farmer = 1/21
pr_librarian = 1 - pr_farmer
prior = [pr_farmer,pr_librarian]#事前確率
posterior = [0.087,1-0.087]#事後確率
plt.bar([0,.7],prior,alpha=0.7,width=0.25,color=colors[0],label="Prior distribution",lw="3",edgecolor="#348ABD")#事前確率
plt.bar([0+0.25,.7+0.25],posterior,alpha=0.7,width=0.25,color=colors[1],label="Posterior distribution",lw="3",edgecolor="#A60628")#事後確率
plt.xticks([0.2,0.95],["Librarian","Farmer"])#司書, 農家
plt.ylabel("Probability")#確率
plt.legend(loc="upper left")
plt.title("Prior and posterior probabilities of Steve's occupation")
plt.show()

今回は P(X|not A)を0.5に設定したがそれを0-1.0に変化させたときにどう変わっていくかを観察する.
つまり「Steveが農家だという情報が分かっているときに, Steveが内向的である確率」P(X|notA)の変化によって, Steveが司書、農家である事後確率はどのように変化するのか?

from IPython.core.pylabtools import figsize
import numpy as np
import matplotlib.pyplot as plt

figsize(12.5,4)

plt.figure()
plt.subplots_adjust(wspace=0.4, hspace=0.6)

colors = ["#348ABD","#A60628"]
pr_farmer = 1/21
pr_librarian = 1 - pr_farmer
prior = [pr_farmer,pr_librarian]#事前確率

a = np.linspace(0,1,11)
print(a)
for i,_ in enumerate(a):
    #Steveが農家であるときに内向的である事後確率
    p_x = 0.95*(1/21) + _ *(20/21)
    pos = (0.95*(1/21))/p_x
    posterior = [pos,1-pos]#事後確率
    print(i)
    plt.subplot(6,2,i+1)
    plt.bar([0,.7],prior,alpha=0.7,width=0.25,color=colors[0],label="Prior distribution",lw="3",edgecolor="#348ABD")#事前確率
    plt.bar([0+0.25,.7+0.25],posterior,alpha=0.7,width=0.25,color=colors[1],label="Posterior distribution",lw="3",edgecolor="#A60628")#事後確率

    plt.title("P(X|notA) = "+str(_))
plt.show()

f:id:umashika5555:20170925172951p:plain
これは P(A|X) = \frac{P(X|A)P(A)}{P(X)} = \frac{0.95*\frac{1}{21}}{0.95*\frac{1}{21}+x*\frac{20}{21}}=\frac{0.95}{0.95+20x}となるため,
x=0のとき「Steveが内向的だという情報を得たときのSteveが司書である」事後確率は1.0となる.
しかし,0.2以降は事後確率があまりかわらないことがわかった.
試しにグラフを書いてみた.

import numpy as np
import matplotlib.pyplot as plt
x =  np.linspace(0,1.,1000)
y = 0.95/(0.95+20*x)
plt.scatter(x,y)
plt.show()

f:id:umashika5555:20170925173414p:plain

【scikit-learn】Ridge回帰【その2】

RidgeCVというモデルがあるらしい.
CVとはCrossVaridationのことで,正則化パラメータαに交差検証を用いる方法である.
これによって最適なαを自動で求めてくれるそう.

reg = linear_model.RidgeCV(alpha= [0.1, 1.0, 10.0])
reg.fit(X_train, y_train)
reg.alpha_

でできる.
デフォルトではleave-out-cross-varidation(LOOCV,一個抜き交差検証)で正則化パラメータの候補集合から1つをテスト事例として抜き,残りをトレーニングデータとして使う.
これを全事例が一回ずつテスト事例となるように検証を繰り返す.

【scikit-learn】Ridge回帰

Ridge回帰は係数 wに対して min_{w}||wX-y||^{2}+\alpha||w||^{2}を満たす wを見つけるというもの.

Ridge回帰による正則化係数 \alphaと係数 wの関係を見てみる
import numpy as np
import matplotlib.pyplot as plt
import sklearn.linear_model

# 10*10の行列を作る 
X = 1. / (np.arange(1, 11) + np.arange(0, 10)[:, np.newaxis])
y = np.ones(10)

Xは

[[ 1  2  3  4  5  6  7  8  9 10]
 [ 2  3  4  5  6  7  8  9 10 11]
 [ 3  4  5  6  7  8  9 10 11 12]
 [ 4  5  6  7  8  9 10 11 12 13]
 [ 5  6  7  8  9 10 11 12 13 14]
 [ 6  7  8  9 10 11 12 13 14 15]
 [ 7  8  9 10 11 12 13 14 15 16]
 [ 8  9 10 11 12 13 14 15 16 17]
 [ 9 10 11 12 13 14 15 16 17 18]
 [10 11 12 13 14 15 16 17 18 19]]

の行列の各要素で1を割った行列なので,左上に行くほど大きく,右下に行くほど小さい値となっている0.5~1.0の行列である.
次に学習を行う.

n_alphas = 200
alphas = np.logspace(-10, -2, n_alphas)#1e-10 ~ 1e-2までを10を底とした指数で表示
coefs = []
for a in alphas:
    ridge = linear_model.Ridge(alpha=a, fit_intercept=False)#モデルの定義
    ridge.fit(X, y)#学習
    coefs.append(ridge.coef_)

ridge.fit(X,y)のところは図にするとおそらくこんな感じ.
f:id:umashika5555:20170925141927p:plain
次にプロットを行う.

ax = plt.gca()
ax.plot(alphas, coefs)
ax.set_xscale('log')
ax.set_xlim(ax.get_xlim()[::-1])  # reverse axis
plt.xlabel('alpha')
plt.ylabel('weights')
plt.title('Ridge coefficients as a function of the regularization')
plt.axis('tight')
plt.show()

出力結果は下のようになった.
f:id:umashika5555:20170925143538p:plain
これはあるαのときにおける式1-10の係数wの大きさを表している.
αが小さいとき(グラフの右に行くほど) min_{w}||wX-y||^{2}+\alpha||w||^{2} = min_{w}||wX-y||^{2}となるのでこれはほぼ線形回帰であり,係数の影響を抑えきれなくなる.
αが大きいほど(グラフの左に行くほど)wは小さくなりたがろうとするので係数の大きさはほぼ0になる.

【補足】

ベクトルの作り方
a = np.arange(0,10)
#[0,1,2,3,4,5,6,7,8,9]
b = np.linspace(0,10,6)
#[0,2,4,6,8,10]
c = np.logspace(2,5,4)
#[10**2,10**3,10**4,10**5]
d = np.ones(10)
#[1,1,1,1,1,1,1,1,1,1]
e = np.zeros(10)
#[0,0,0,0,0,0,0,0,0,0]

【参考】
http://scikit-learn.org/stable/auto_examples/linear_model/plot_ridge_path.html#sphx-glr-auto-examples-linear-model-plot-ridge-path-py

【numpy】行列操作【その1】

numpyの行列操作で面白い行列の作り方を知ったので書く.

a = np.arange(1,11)
b = np.arange(0,10)
#b: [0,1,2,3,4,5,6,7,8,9]
c = b[:,np.newaxis]
#c: 
#[[0]
# [1]
# [2]
# [3]
# [4]
# [5]
# [6]
# [7]
# [8]
# [9]]
d = a + c
#[[ 1  2  3  4  5  6  7  8  9 10]
# [ 2  3  4  5  6  7  8  9 10 11]
# [ 3  4  5  6  7  8  9 10 11 12]
# [ 4  5  6  7  8  9 10 11 12 13]
# [ 5  6  7  8  9 10 11 12 13 14]
# [ 6  7  8  9 10 11 12 13 14 15]
# [ 7  8  9 10 11 12 13 14 15 16]
# [ 8  9 10 11 12 13 14 15 16 17]
# [ 9 10 11 12 13 14 15 16 17 18]
# [10 11 12 13 14 15 16 17 18 19]]
X = 1./ d
# 1.をdの要素で割った形

【参考】
http://scikit-learn.org/stable/auto_examples/linear_model/plot_ridge_path.html#sphx-glr-auto-examples-linear-model-plot-ridge-path-py

matplotlib.pyplotで一枚の枠に図を複数埋め込む方法

下のように宣言すればよい

plt.figure()
plt.subplot(y,x,i)

これは文字列でもできるらしく,scikit-learnの公式documentではこちらの方が使用されていたため,文字列のほうが主流なのか?

plt.subplot(211)
plt.subplot(212)
正規分布を作る

numpy.random.randn()は正規分布に従って乱数が生成されるらしい.
今回は
x = numpy.random.randn(発生させる個数)
y = 正規関数(x)
正規分布の曲線の発生過程を作った.

import matplotlib.pyplot as plt
import numpy as np

plt.figure()
h_length = 4
w_length = 6
for h in range(h_length):
    for w in range(w_length):
        i = h*w_length+w+1
        fig_loc = str(h_length) + str(w_length) + str(i)
        plt.subplot(h_length,w_length,i)
        #グラフプロット
        #xの個数
        num = 10*i
        x = np.random.randn(num)
        y = np.exp(-x**2/2)/np.sqrt(2*np.pi)
        plt.title(num)
        plt.xlim([-3,3])
        plt.ylim([-0.1,0.55])
        plt.scatter(x,y,c="green")
plt.show()

f:id:umashika5555:20170924235757p:plain

【scikit-learn】線形回帰

研究で識別/分類問題について実装することになるだろうので,scikit-learnの公式ドキュメント(http://scikit-learn.org/stable/index.html)のうち
classificationを勉強していこうと思う.
今回は線形回帰について(なぜ回帰?)

scikit-learnは

#modelの定義
regression = sklearn.linear_model.LinearRegression()

#学習
#定義したモデル.fit(トレーニングデータ, トレーニングデータのラベル)
regression.fit(X_train,y_train)

#予測
#予測ラベル = 定義したモデル.predict(テストデータ)
y_pred = regression.predict(X_test)

という感じでやるらしい.
線形回帰はMSEを最小化するという方法である.

糖尿病患者

diabetesデータセット
age, sex, bmi, map, tc, ldl, hdl, tch, ltg, gluという10個の属性が入っている.
目的変数は1年後の疾患の進行状況らしい.
これらを線形回帰によって当てはめてみる.

import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets
from sklearn import linear_model
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score

# load the diabete dataset
diabetes = datasets.load_diabetes()
categories = {0:"age",1:"sex",2:"bmi",3:"map",4:"tc",5:"ldl",6:"hdl",7:"tch",8:"ltg",9:"glu"}

#全体の図を定義
plt.figure()

for key,val in categories.items():
    print(key)
  #データをある属性のみに絞る
    diabetes_X = diabetes.data[:, np.newaxis, key]
    #データをトレーニングとテストに分割
    diabetes_X_train = diabetes_X[:-20]
    diabetes_X_test  = diabetes_X[-20:]
    #目的変数をトレーニングとテストに分割
    diabetes_y_train = diabetes.target[:-20]
    diabetes_y_test  = diabetes.target[-20:]
    #線形回帰モデルの宣言
    regr = linear_model.LinearRegression()
    #学習
    regr.fit(diabetes_X_train,diabetes_y_train)
    #予測
    diabetes_y_pred = regr.predict(diabetes_X_test)
    #係数
    print("Cofficients:\n",regr.coef_)
    #MSE計算
    print("mean squared error:%.2f"%(mean_squared_error(diabetes_y_test,diabetes_y_pred)))
    #分散
    print("Variance score:%.2f"%r2_score(diabetes_y_test,diabetes_y_pred))
    #subplot
    fig_loc = "25" + str(key)
    plt.subplot(fig_loc)
    #plot outputs
    plt.scatter(diabetes_X_test,diabetes_y_test,color="black")
    plt.plot(diabetes_X_test,diabetes_y_pred,color="blue",linewidth=3)
    #plt.xticks(())
    #plt.yticks(())
    plt.title(categories[key])
plt.show()

結果はこのようになった.
sexのような2値に対しては当然ながら線形回帰は使えない.(そもそも使う意味がない.)
bmiやhdl, ltgあたりは割りと線形回帰で良い結果が出た方だと思う.
f:id:umashika5555:20170924230252p:plain

識別器に意図的に誤分類させる研究

識別器に意図的に誤分類を誘発させる研究の記事を読んだ.
今のプロジェクトが終わったらやってみたいのでメモ.

TeX余白の設定

箇条書における余白の設定

http://sumi-nyasu.hatenablog.com/entry/2015/06/27/232345

表や図における余白の設定

http://d.hatena.ne.jp/pyopyopyo/20120416/p1