Pythonで積分計算

Python

区分求積法より
\sum_{k=0}^{n-1}\Delta xf(\frac{k}{n}) = \frac{1}{n}\sum_{k=0}^{n-1} f(\frac{k}{n})
と
\sum_{k=0}^{n}\Delta xf(\frac{k}{n}) = \frac{1}{n}\sum_{k=0}^{n-1} f(\frac{k}{n})
をはさんだものが積分の値となるので分割数nを大きくすればその値に近似できる.
今回は y = (1-x)^{2}xについて計算する.
まずは区分求積法をnumpyで実装.

import numpy as np
import matplotlib.pyplot as plt
n = 1000
x = np.linspace(0,1,n)
y = x*(1-x)**2
ans = np.sum(y)/n
print(ans)#0.0832499165832

科学計算用のライブラリが用意されている.

import nunmpy as np
import scipy import integrate

def compute(x):
    return x*(1-x)**2

ans = integrate.quad(compute,0,1)# ans = integrate.quad(lambda x:x*(1-x)**2,0,1)
print(ans)#(0.08333333333333334, 9.251858538542972e-16)

こちらは誤差まで議論されている.

はてなブログのデザインを一部いじる

まずオリジナルのデザインはこちら
f:id:umashika5555:20170927015532p:plain
今回はブラウザがGoogleChromeでやることを想定する.

1.F12キーを押すとブラウザの右か下にバーが現れるので, その左上マウスのマークをクリック
f:id:umashika5555:20170927015629p:plain

2.オリジナルのデザインで変更したい部分にマウスを当てクリックする

3.該当箇所のCSSが表示される
今回は

.entry-content pre.code {
    color: #333;
    border: none;
    padding: 20px;
    background: #f6f6f2;
    font-size: 85%;
    margin: 0.8em 0;
}

のbackgroundを#f6f6f2から#CCCに変更する.
このCSSを全てコピーする.

4.はてなのメニューの「デザイン」から「デザインCSS」をクリック
そこに先ほどのコピーしたCSSを貼り付ける.
改変したい箇所#f6f6f2を#CCCに変更する.

5.デザイン変更完了
f:id:umashika5555:20170927020400p:plain

【Ubuntu】MeCabインストール

Python

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

【Windows】TeXworksにforestを入れる

Windows TeX

dirtreeパッケージよりもforestのほうが綺麗にディレクトリツリーを表示できたのでメモ.

インストールの仕方

forestパッケージのダウンロード

https://www.ctan.org/pkg/forest
https://github.com/sasozivanovic/forest

inlinedefパッケージのダウンロード

forestだけで宣言したらinlinedef.styがないとエラーが出たので依存関係がある
https://www.ctan.org/tex-archive/macros/latex/contrib/inlinedef

ビルドする

各ディレクトリに行って

platex *.ins

TeXworksのディレクトリに入れる

自分の環境だとC:\w32tex\share\texmf-dist\tex\latex\0_mypackage\

TeXの環境を更新

$mktexlsr

サンプル

\usepackage{forest}
%forest definition-------------------------------------------
\definecolor{folderbg}{RGB}{124,166,198}
\definecolor{folderborder}{RGB}{110,144,169}

\def\Size{4pt}
\tikzset{
  folder/.pic={
    \filldraw[draw=folderborder,top color=folderbg!50,bottom color=folderbg]
      (-1.05*\Size,0.2\Size+5pt) rectangle ++(.75*\Size,-0.2\Size-5pt);  
    \filldraw[draw=folderborder,top color=folderbg!50,bottom color=folderbg]
      (-1.15*\Size,-\Size) rectangle (1.15*\Size,\Size);
  }
}
%-------------------------------------------
\begin{document}
\begin{figure}

 \begin{forest}
      for tree={
        font=\ttfamily,
        grow'=0,
        child anchor=west,
        parent anchor=south,
        anchor=west,
        calign=first,
        inner xsep=12pt,
        edge path={
          \noexpand\path [draw, \forestoption{edge}]
          (!u.south west) +(6.5pt,0) |- (.child anchor) pic {folder} \forestoption{edge label};
        },
        % style for your file node 
        file/.style={edge path={\noexpand\path [draw, \forestoption{edge}]
          (!u.south west) +(6.5pt,0) |- (.child anchor) \forestoption{edge label};},
          inner xsep=2pt,font=\small\ttfamily
                     },
        before typesetting nodes={
          if n=1
            {insert before={[,phantom]}}
            {}
        },
        fit=band,
        before computing xy={l=15pt},
      }  
    [Project
    	[Sample
    		[sample.py,file]
    		[sample.java,file]
	]
	[img
		[a.png,file]
		[b.png,file]
	]
    ]
\end{forest}
\centering\caption{プロジェクトの構成}\label{Project}
\end{figure}
\end{document}

f:id:umashika5555:20170913220033p:plain

WidowsのTeXworksでstyファイルを読み込む

TeX

今回はdirtreeというパッケージを入れた.
忘れないようにメモ

TeXの\usepackage{hoge}という宣言ではhoge.styというファイルを読み込んでいる.
最も簡単な方法はそのtexファイルと同じディレクトリに入れることである.
しかしそのプロジェクトのみでしか使えないため不便である.
そこでPATHの設定のようにいつでも使えるように設定した.


1. まずdirtreeファイル一覧をダウンロードする.
https://ctan.org/pkg/dirtree


2. insファイルからstyファイルを生成する

$ platex dirtree.ins


3. TeXworksのフォルダに入れる
自分の場合だとC:\w32tex\share\texmf-dist\tex\latex\0_mypackage\dirtreeというフォルダにファイル一覧を突っ込んだ.
先に読み込まれるほど優先度が高いそうなのでファイル名は0などから始めると良いらしい.

4. 環境を更新

$ mktexlsr


5. 確認

\usepackage{dirtree}

でエラーが出ないか確認

【参考】
http://oku.edu.mie-u.ac.jp/tex/mod/forum/discuss.php?d=454
https://ameblo.jp/masayan24/entry-10360005851.html
http://ocg.aori.u-tokyo.ac.jp/member/daigo/comp/memo/?val=all&typ=all&nbr=0000000170