2016-05-04

RをJupyter notebookで動かす

R Jupyter

Jupyter notebookでRを動かしたいと思ったのですが、どうにも公式のGitHub 通りにやっても上手くいかない。

ずーっと詰まっていたのですが、結局、これにしたがってやるとすんなりうまく行った。
www.continuum.io

Jupyterは、コードを管理したり、共有するのがすごく得意と聞いたので、素人ながら、試しにコードを共有してみます。

最近、ヨットの関東大会の成績表を勝手にWebでつくったのですが、その際のページを共有します。

https://web.sfc.keio.ac.jp/~t13804kf/haruin2016/

2015-10-11

ブンセキマラソンその１　「要約統計量・度数集計・ヒストグラム」

始まりましたブンセキマラソン(仮)初回は要約統計量・度数集計・ヒストグラムの３点盛り！

説明のためにまずはヒストグラムから解説します。

ヒストグラム

f:id:fujit33:20151003214102j:plain
ヒストグラムとはこのような棒グラフっぽいものでデータの分布を視覚的に捉えるものです。でも棒グラフと違うのは、値を区切ってまとめていること！
たとえば、このヒストグラムでは130と140の間に3人居ることを表しています。ここでは「130cm台」としてまとめてしまっていますが、本当は133,136,139というデータが中にあります。このように、データをまとめることで全体の傾向が見やすくなっているのがヒストグラム。
ヒストグラムは棒の間に隙間を作らずに描きます。

要約統計量とは

データを要約した値が要約統計量。例を挙げたほうが早い。

平均(mean)

データの合計をデータの個数で割ったもの。とってもよく使う。
${ \displaystyle 平均 \mu = \frac{1}{N}(x_1+x_2+x_3\cdots) = \frac{1}{N}\sum_{i=1}^{N} x_i }$

中央値(median)

得られたデータを小さい順に並べて、ちょうど真ん中の順番になった値。ただし、データが偶数の場合、真ん中の２つの値の平均値。

例：
「0,1,2,3,4,5,6,7,8,9,10」というデータの中央値は「5」

最頻値(mode)

データの中で回数が最も多く現れる値。複数あることもある。

例：
「0,1,2,3,3,3,4,4,5,5,6」というデータの最頻値は「3」

分散・標準偏差

分散・標準偏差はデータの「ばらつき」を示す値。標準偏差は分散の平方根（ルート）なので、それぞれが示す意味は同じ。

分散はデータの値ひとつひとつの「平均との差」を二乗し（マイナスの差もプラスにするため）てすべて合計し、それをデータの個数で割った（データの個数に影響されないため）もの。
式で表すと
${ \displaystyle 分散V = \frac{1}{N}{(x_1-\overline{x})^2+(x_2-\overline{x})^2+(x_3-\overline{x})^2\cdots} = \frac{1}{N}\sum_{i=1}^{N} (x_i-\overline{x})^2 }$
${ \displaystyle 標準偏差σ = \sqrt{V} }$
分散や標準偏差が大きいほど、平均から離れたところにデータが沢山あるということ。

例えば、同じ人たちが３つのテストを受験したらどれも平均は60点だった。しかし、標準誤差が違った場合、このようなヒストグラムになります。
f:id:fujit33:20151005151428j:plain
赤いテストは標準偏差10、緑は15、青は20
標準偏差が小さいほど、データがばらつかず、標準偏差が大きいほどデータがばらけることがよく分かると思います。

平均信じるな！

平均値はとてもわかりやすく便利な値ですが、むやみに平均値を鵜呑みにするとデータを間違って解釈する恐れがあります。実際、世の中のデータは平均値で表せないデータばかりです。
例えば、次のグラフはTwitterにおける、一日あたりの平均ツイート数のヒストグラムです。
f:id:fujit33:20151011202855j:plain
データがかなり左に偏っているが、右にも長く伸びていることがわかります。つまり、ほとんどの人は１日あたりツイート数が少ないが、少数の人々は日々大量にツイートを投稿するヘビーユーザーであるとわかります。このようなデータで中央値と平均値を計算してみるとどうなるでしょうか。
f:id:fujit33:20151011202845j:plain
中央値と平均値にかなり大きな差がでてしまうことが分かります。このように、大きく偏っているデータでは、平均値よりも中央値を用いたほうが良いことがあることを覚えておかなければいけません。（注意点として、このデータは機械的に自動投稿するユーザーも含んでいるので、感覚よりも多くなっています。）

＜参考文献＞
http://www.agri.tohoku.ac.jp/iden/toukei2.html

2015-10-02

分析手法マラソン始めます

その他

ず～～～っと更新せずに秋が来ました

最近、自分に何が足りていないのか考えることが良く有ります。それはきっと就活の時期がひしひしと近づいてきているから。
そして決心しました。マラソンをしようと。

これまで大学でデータ分析についていろいろ分析して、テーマを決めて研究もして、データ分析コンテストに出たりもして、いろいろやって来ました。その結果今の自分を分析してみると、自分が自信を持ってアピールできるポイントがひとつもないことに気づいたのです。分析もちょっと手法を知ってて、ちょっとRを使えてちょっと応用できるだけ。本質的に分析手法を理解して使いこなせて無いし、これじゃあ大学の授業の実習となんら変わらない。

理論を学び ➔ 応用するもっともっとこの経験を積まないと１流にはなれない。院卒には勝てない。そう思ったのです。

ということで、一通りの分析手法について、ブログ上で理論を解説し、何かしらの分析をする分析手法マラソンを一人で開催します！

記念すべき第一回大会（第二回はやらない、絶対に）は『データサイエンス超入門ビジネスで役立つ「統計学」の本当の活かし方』という本の付録にあった「構造化データサイエンスモデル」に書かれている分析手法を順番に、週一回のペースで解説、分析していきます。

f:id:fujit33:20151002234733p:plain ]

初回は丁寧に要約統計量からやり、ヒストグラムまでやります。

このマラソンを完走したとき、自分が大きな成長を遂げていることを信じています。

分析手法マラソンって名前カッコ悪いからもっといい名前無いかな。

ツール・ド・データブンセキ…………ありだな

2015-05-14

暖かくなってきたので、GitHub始めました

その他

最近、暖かい日、暑い日が多くなってきました。というわけでGitHubを始めました。

正直、人に見せられるコードは書いていないのですが、やっぱり誰かに見せるつもりで書かないと成長しないと思い、始めました。（真の理由）

<a href="https://github.com/fujit33">fujit33</a>github.com

じゃんじゃんダメ出しお願いします。

2015-05-14

ナイーブベイズ分類器でTwitterの男女判別！

R Twitter 機械学習

Twitterの分析をする上で避ける事ができないのが男女判別。Twitter Analiticsではすでに自分のフォローワーの男女比が表示されます。
Twitter者はフォロー関係などを使って判別しているらしいのですが、今回はつぶやかれているテキストを使って、「ナイーブベイズ」という手法で分類してみます。

データについて

そもそも今回は勉強中のPythonでツイート収集、分析を行おうと思ったのですが、途中でどうしても解決できないエラーが有ったので一旦諦めました。
データは性別を目視で確認できた男62ユーザー,女93ユーザー（両方100人くらい集めたはずが途中で闇に消えました……）からそれぞれの最新1000ツイートを取得しテキストファイルとして保存してあります。

ナイーブベイズとは

ナイーブベイズとは「単純ベイズ分類器」とも呼ばれ、スパムメールフィルターなどに用いられています。理論はシンプルながらも、高い精度を出すことができます。（数式が出てきますが、頑張ればきっと理解できるように書きました？）

ナイーブベイズの理論はベイズの定理を応用して次のように表せます。

${ \displaystyle P(cat|doc) = \frac{P(doc|cat)P(cat)}{p(doc)} \propto P(doc|cat)P(cat) }$

$P(cat|doc)$ は文章docが与えられた時にカテゴリcatに属する確率（これを知りたい！）
$P(doc|cat)$ はカテゴリcataに属するときに、文章docである確率（わかりづらいのであとで解説）
$P(cat)$ はカテゴリcatである確率（catの文章数 / 全文章数)
$P(doc)$ は文章docである確率

$\propto$ は右に比例するという意味です。
今回は確率そのものを知りたいわけではなく、最も確率が高いカテゴリを知りたいので、確率に比例する $P(doc|cat)P(cat)$ を最大化すればいいわけです！

具体的に考える

まず、 $P(cat)$ は、男(女)の文書数 / 全体の文書数で簡単に出せます。

次に $P(doc|cat)$ ですがけっこうややこしいです。

ナイーブベイズでは、まず文章を単語に切り分け、単語の集合を「文書」として扱います。集合なので、文章の中で出てくる順番は関係なく、このような表現をbag-of-wardsと呼ばれます。
そしてそれぞれの単語の出現回数を文章ごとに調べてこのような行列にします。
f:id:fujit33:20150514012713p:plain:w500
この行列を使うと、単語Aの文書T内での出現回数 / 単語Aの総出現回数　を調べることができます。
それをカテゴリに属する文章の、全ての単語に関してしらべ、掛け合わせたものが $P(doc|cat)$ です！！！！
数式にすると

$P(doc|cat) = P(word_i|cat) \cdots = \prod_i P(word_i|cat)$

となります。 $\prod$ は $\sum$ の掛け算バージョンで、iが1ずつ増えていき、単語を順番に変えながら確率を出し、それらを全て掛け合わせると言う意味です。

ついに求める数式が！？

これらより、最終的に求めるカテゴリmapは

$cat_{map} = arg max_{cat} P(doc|cat) = arg max_{cat} P(cat)\prod_i P(word_i|cat)$

となります。 $arg max_{cat}$ とはcatについて最大化する値を返すということです。しかし問題が有り、文章内の単語数は膨大なため、 $P(word_i|cat)$ をすべてかけ合わせるとひっじょーーーーに小さい数字になり、計算がアンダーフローを起こしてしまう可能性があります。（あるそうです）アンダーフローとは、簡単にいえば数字が小さくなりすぎて、そのプログラミング言語じゃ処理できないよってことです。

だからといってここまできて諦めるわけにはいかないので、対数をとって掛け算ではなく足し算の計算に変えてあげます。「対数の計算なんて忘れた」って人はググッて下さい。

すると式はこう変換できます。

$cat_{map} = arg max_{cat} P(doc|cat) = arg max_{cat} \big(log(P(cat)) + \sum_i log(P(word_i|cat))\big)$

ラプラススムージング

実はこの式のままでは欠点があって、 $log(P(word_i|cat))$ は単語が一回も出現しない場合だと0になってしまい、log 0は計算できません。それを回避するために、すべての単語に出現回数を１足してあげるという処理をします。この処理をラプラススムージングと呼びます。

Rで実装

ここまでをようやくRで実装するのですが、今回はあらびきさん(@a_bicky)のこちらのコードをお借りしました。少し解説を加えています。10行でナイーブベイズを実装してしまうあらびきさんはすごいです……

library(RMeCab)

d <- t(docMatrix2("usertimelines"))

ncol(d)
myNaiveBayes <- function(x, y) {
    # 文章のカテゴリ
    lev <- levels(y) #1
    # それぞれのカテゴリでの単語の出現数
    ctf <- sapply(lev, function(label) colSums(x[y == label,])) #2
    # ラプラススムージングを行ったそれぞれのカテゴリでの単語の出現確率
    ctp <- t(t(ctf + 1) / (colSums(ctf) + nrow(ctf))) #3
    # それぞれのカテゴリの文章数
    nc <- table(y, dnn = NULL) #4
    # カテゴリの文書数の割合
    cp <- nc / sum(nc) #5
    structure(list(lev = lev, cp = cp, ctp = ctp), class = "myNaiveBayes") #6
}

predict.myNaiveBayes <- function(model, x) {
    # モデルにそれぞれの単語の出現確率の対数をすべて足す
    prob <- apply(x, 1, function(x) colSums(log(model$ctp) * x)) #7
    # 学習データの文書量で重みを変える
    prob <- prob + log(as.numeric(model$cp)) #8
    # 確率が高い方のインデックスを返す
    level <- apply(prob, 2, which.max) #9
    # 確率が高いほうのラベルを返す
    model$lev[level] #10!!
}


y <- factor(sub("^([a-z]*?)_.*", "\\1", rownames(d), perl = TRUE)) # データにラベル付け
table(y)
train.index <- c(1:42, 63:135) # 学習データにする文書の選択
model <- myNaiveBayes(d[train.index,], y[train.index]) # 学習
test <- y[-train.index] # テストデータの回答ラベル
result <- predict(model, d[-train.index,]) # テスト

# 結果発表(ドキドキ）
table(test, result) 
mean(test == result)

結果

> table(test, result)
       result
test    man woman
  man    18     2
  woman   2    18
> mean(test == result)
[1] 0.9

正解率は90%でした。

まとめ

Twitterは特殊な文章であるので、全然分類できないことを覚悟していたので案外すんなりいってうれしいですが、肩透かし感もあります。90%分類できたのですが、まだ工夫によって正解率はあげられるはずです。例えば、絵文字や顔文字など、Twitter特有のよく使われる文字に注目すると効果があると思います。
また、今回は男女の２分類でしたが、カテゴリを増やした分類もやってみたいと思います。どうやら、多項モデルの方が精度が上がるという論文もあるようです。でも今はいろいろな手法を試したいから後回しかな……

どうでもいいつぶやき
最近、名著『パターン認識と機械学習』(PRML)で機械学習の勉強をしています。難しい……

(参考リスト)
単純ベイズ分類器 - Wikipedia
第3回ベイジアンフィルタを実装してみよう：機械学習はじめよう｜gihyo.jp … 技術評論社
 ナイーブベイズを用いたテキスト分類 - 人工知能に関する断創録

2015-05-08

MacのSublime Text 3 で行頭、行末に移動するキーバインドの設定

Mac Sublime Text

普段エディタで、Sublime Text 3を使っているのですが、Linux風のctrl + e(a)で行末（行頭）に移動するキーバインドを使えるようにするのに苦労したのでメモ。

そもそもデフォルトではctrl + e(a) をすると文章の最後（最初）になってしまう。ネットを探しまわったのですが、同じ症状の人は見当たらなかったので、もしかしたら私だけかも、、、

キーバインドを設定するには
Sublime Text　→　Preferences →　Key Binding User
を編集します。

SublimeText3で行頭・行末に移動するキーボードショートカット | DevAchieve
ここを参考に

  // 行末に移動
  { "keys": ["ctrl+e"], "command": "move_to", "args": { "to": "hardeol" } },
  // 行頭に移動
  { "keys": ["ctrl+a"], "command": "move_to", "args": { "to": "hardbol" } },

を[　]の間に書き込みます。

が...

何も変わりませんでした。
ここでずっと止まって居たのですが、ふとeやaじゃなくて別のアルファベットにすればどうなるかと試してみると、なんとできたのです。うーん分からない。でも、それでできるならいいやって感じです。

現状はこんな感じ。

  // 行末に移動
  { "keys": ["ctrl+m"], "command": "move_to", "args": { "to": "hardeol" } },
  // 行頭に移動
  { "keys": ["ctrl+c"], "command": "move_to", "args": { "to": "hardbol" } },

2015-04-12

Rで目にやさしい色を使う

rでプロットするとき、col="red"とかcol="blue"とかやるが、なんとも目に痛い色で嫌だ。

お気に入りの自分用パレットをつくろうかと思ったが、RColorBrewerパッケージが解決してくれた。

library(RColorBrewer)
cols <- brewer.pal(8,"Set1") # 引数1:何色使うか　引数2:パレットの名前

あとは色を指定するとき、cols[1]やcols[8]とするだけで、8色の目にやさしい色が使える。ちなみに.Rprofileに書いておくと便利。

例えばこんな感じで使う

plot(iris$Sepal.Length,iris$Sepal.Width,pch=16,col="red")
plot(iris$Sepal.Length,iris$Sepal.Width,pch=16,col=cols[1])
plot(iris$Sepal.Length,iris$Sepal.Width,pch=16,col="blue")
plot(iris$Sepal.Length,iris$Sepal.Width,pch=16,col=cols[2])
plot(iris$Sepal.Length,iris$Sepal.Width,pch=16,col="green")
plot(iris$Sepal.Length,iris$Sepal.Width,pch=16,col=cols[3])

f:id:fujit33:20150412210611p:plain:w300 f:id:fujit33:20150412210618p:plain:w300
f:id:fujit33:20150412210612p:plain:w300 f:id:fujit33:20150412210615p:plain:w300
f:id:fujit33:20150412210614p:plain:w300 f:id:fujit33:20150412210616p:plain:w300