JISキーボードの記号の並び

　アスキーコード表を眺めていて見たことある並びがあるなと思ったら、JISキーボードの並びと同じだった。

　JISキーボードでは、!"#$%&'()の記号が123456789のキーに配置されている。
!"#$%&'()のアスキーコードは21,22,23,24,25,26,27,28,29（16進数表記）であり、連番である。

　つまりJISキーボードを見ると、これら9つの記号のアスキーコードの大小がすぐに分かる。いや、それだけではない。アスキーコードを16で割った余りさえも分かってしまうのだ（キーの数字を見ればよい）。

というどうでもいい話に気づいてしまった2016年の正月早朝....

2015年を画像一枚で振り返る

 Pythonのwordcloudというライブラリを使って今年一年を振り返ってみました。

今年の振り返り

今年書いたブログの内容をもとに作った画像が上です。
何だ、このリア充のリの字も感じさせない結果は....
来年は「女」「酒」「車」とかが入ってくるよう精進したいです。

ソースコード

# -*- coding: utf-8 -*-

import time
import os
import re

import urllib3
from bs4 import BeautifulSoup
import MeCab
import matplotlib.pyplot as plt
from wordcloud import WordCloud

class YearWordVisualizer():
    def __init__(self, blog_url, year):
        self.blog_url = blog_url
        self.year = year

    def run(self):
        self.crawl_articles()
        self.preprocess()
        self.visualize()
    
    def crawl_articles(self):
        http = urllib3.PoolManager()
        for month in range(1, 13):
            time.sleep(1)
            url = '{blog}/{year}/{month:02d}'.format(blog=self.blog_url, year=self.year, month=month)
            r = http.request('GET', url)
            open('html/{month:02d}.html'.format(month=month), 'w').write(r.data)
    
    def preprocess(self):
        cnt = 0
        for f_name in os.listdir('html'):
            with open('html/'+f_name, 'r') as f:
                html = f.read()
                soup = BeautifulSoup(html, 'lxml')
                for pre in soup.find_all(['pre', 'img']):
                    pre.decompose()
                for post in soup.find_all('div', class_='post-body'):
                    text = post.getText()
                    nouns = self.get_noun_list(text.encode('utf-8'))
                    doc = " ".join(self.remove_stop_words(nouns))
                    cnt += 1
                    output_file = "{dir}/{cnt}.txt".format(dir='text', cnt=cnt)
                    open(output_file, 'w').write(doc)
                    
    def get_noun_list(self, text):
        tagger = MeCab.Tagger('-Ochasen')
        node = tagger.parseToNode(text)
        res = []
        while node:
            if node.feature.split(',')[0] == '名詞':
                res.append(node.surface)
            node = node.next
        return res
    
    def remove_stop_words(self, words):
        words = [w for w in words if not w.isdigit()]
        words = [w for w in words if not (w.isalpha and len(w) == 1)]
        pattern = re.compile(r'^[+-/%\^\[\]\(\)\{\}「」（）［］。、]+$')
        words = [w for w in words if not re.match(pattern, w) ]
        words = [w for w in words if w not in ['こと', 'これ', 'それ', 'もの', 'の', 'よう', 'とき', '以下', '場合', '結果']]
        return words
    
    def visualize(self):
        text = ''
        for f_name in os.listdir('text'):
            text += ' '
            text += open('text/'+f_name, 'r').read()
        font_path = '/System/Library/Fonts/ヒラギノ明朝 ProN W3.otf'
        wordcloud = WordCloud(font_path=font_path, width=900, height=500).generate(text.decode('utf-8'))
        plt.imshow(wordcloud)
        plt.axis("off")
        plt.show()

if __name__ == '__main__':
    visualizer = YearWordVisualizer('http://techtipshoge.blogspot.jp', 2015)
    visualizer.run()

LDAで日経225企業をクラスタリングする

　最近LDAという言葉をよく聞くので、LDAを使って何かをやってみることにした。
とりあえず使ってみることを重視するため、深い/細かい話には立ち入らない。

LDAとは？

• Latent Dirichlet Allocationの略。
• トピックモデルの一種。
• ドキュメントは複数のトピックから成り立っていると仮定。
• 入力=ドキュメント集合、出力=各ドキュメントのトピック分布、各トピックの語句分布

やること

日経225の企業をwikipediaの文章を元にクラスタリングする。

事前準備として、225銘柄一覧に記載されている企業のwikiページのHTMLを取得しておく。HTMLはhtml/配下に企業名.htmlというファイル名で格納しておく。

$ ls html | head -n 10
ANAホールディングス.html
DOWAホールディングス.html
IHI.html
J.フロント リテイリング.html
JFEホールディングス.html
JXホールディングス.html
KDDI.html
MS&ADインシュアランスグループホールディングス.html
NTN.html
NTTデータ.html

前処理

html/配下に格納されたHTMLを前処理して、text/配下に格納する。

行った前処理は以下のとおり。

• Javascript/CSSの記述削除
• HTMLタグの除去
• MeCabで分かち書き

# -*- coding: utf-8 -*-

import os
from bs4 import BeautifulSoup
import MeCab

for f_name in os.listdir('html'):
    f_path = 'html/' + f_name
    with open(f_path, 'r') as f:
        data = f.read()
        soup = BeautifulSoup(data, 'lxml')
        for script in soup.find_all('script'): script.decompose()
        for script in soup.find_all('style'): script.decompose()
        text = soup.getText()
        text = [line for line in text.splitlines() if line]
        text = '\n'.join(text)

        tagger = MeCab.Tagger('-Owakati')
        wakati_text = tagger.parse(text.encode('utf-8'))    

        t_path = 'text/{name}.txt'.format(name = os.path.splitext(f_name)[0])
        open(t_path, 'w').write(wakati_text)

LDAモデルの学習
前処理したテキストを使ってLDAモデルを学習する。

• Dictionaryの作成（tokenのid化、document ごとのtokenの頻度計算など）
• 汎用的な語句（60%以上のドキュメントで使われる語句）の除去
• bag of words形式の行列作成
• LDAモデルの構築
• ドキュメント-トピック配分の行列を作成し保存

import logging
import os
import pickle
import gensim

logging.basicConfig(format='%(asctime)s : %(levelname)s : %(message)s', level=logging.INFO)
TOPIC_NUM = 30

texts = [file('text/' + f_name, 'r').read().split() for f_name in os.listdir('text')]
dictionary = gensim.corpora.Dictionary(texts)
dictionary.filter_extremes(no_below=1, no_above=0.6, keep_n=None)
corpus = [dictionary.doc2bow(text) for text in texts]
lda = gensim.models.ldamodel.LdaModel(corpus, num_topics=TOPIC_NUM, id2word = dictionary, passes=20)

topics = []
for c in corpus:
    topic = [0] * TOPIC_NUM
    for (tpc, prob) in lda.get_document_topics(c):
        topic[tpc] = prob
    topics.append(topic)

pickle.dump(topics, open("topics.p", "w"))

クラスタリング
トピック配分を特徴量としてK−meansでクラスタリング。

• 特徴量のL2ノルム正規化
• K-meansの学習
• ドキュメントのクラスタリング

# -*- coding: utf-8 -*-

import pickle
import os
from sklearn.preprocessing import Normalizer
from sklearn.cluster import KMeans

N_CLUSTER = 20

xs = pickle.load(open('topics.p', 'r'))
ys = [os.path.splitext(f_name)[0] for f_name in os.listdir('text')]

xs = Normalizer().fit_transform(xs)
kmeans = KMeans(n_clusters=N_CLUSTER, init='k-means++', max_iter=100, n_init=1)
kmeans.fit(xs)
clusters = kmeans.predict(xs)

for i in range(N_CLUSTER):
    text = 'cluster {clr}:'.format(clr = i)
    companies = [y for k, y in enumerate(ys) if clusters[k] == i]
    print text + " ".join(companies)

結果
それっぽく分かれているところを色つけしてみた。

cluster 0: SUMCO りそなホールディングス オークマ クラレ コニカミノルタ サッポロホールディングス トクヤマ ユニチカ 三井化学 三井造船 三井金属鉱業 三菱マテリアル 双日 商船三井 大林組 宇部興産 帝人 日新製鋼 日本製紙 日本製鋼所 日東電工 日立造船 昭和電工 東レ 東宝 => 三井グループ

cluster 1: JXホールディングス MS&ADインシュアランスグループホールディングス アステラス製薬 アドバンテスト 三井住友トラスト・ホールディングス 信越化学工業 国際石油開発帝石 大和証券グループ本社 太平洋セメント 損保ジャパン日本興亜ホールディングス 日本水産 日産化学工業

cluster 2: KDDI NTTドコモ アマダホールディングス ディー・エヌ・エー デンソー ブリヂストン 京王電鉄 宝ホールディングス 富士通 東海カーボン => 携帯電話関連

cluster 3: オリンパス キッコーマン セブン&アイ・ホールディングス マルハニチロ 三井物産 明治ホールディングス 清水建設 王子ホールディングス

cluster 4: IHI J.フロント リテイリング ふくおかフィナンシャルグループ ニチレイ 三菱UFJフィナンシャル・グループ 三菱商事 三越伊勢丹ホールディングス 千代田化工建設 第一生命保険

cluster 5: NTTデータ TDK いすゞ自動車 アルプス電気 クレディセゾン コナミホールディングス ソフトバンク テルモ トヨタ自動車 トレンドマイクロ マツダ ヤマトホールディングス 三菱地所 大平洋金属 大成建設 太陽誘電 富士重工業 小田急電鉄 新日鐵住金 日本碍子 日本通運 日野自動車 昭和シェル石油 東武鉄道 東海旅客鉄道 松井証券 沖電気工業 積水ハウス 花王 長谷工コーポレーション => 自動車関連、鉄道関連

cluster 6: ミネベア 中外製薬 味の素 大日本住友製薬 富士電機 日本曹達 日本精工 日清製粉グループ本社 横河電機 武田薬品工業

cluster 7: ジェイテクト スカパーJSATホールディングス 京成電鉄 古河機械金属 古河電気工業 日本軽金属ホールディングス 東日本旅客鉄道 東洋製罐グループホールディングス 西日本旅客鉄道

cluster 8: TOTO クボタ ヤマハ ユニーグループ・ホールディングス 中部電力 日清紡ホールディングス 東京電力 横浜ゴム 関西電力

cluster 9: カシオ計算機 キヤノン ミツミ電機 凸版印刷 富士フイルムホールディングス 川崎汽船 日本化薬

cluster 10: SCREENホールディングス みずほフィナンシャルグループ コムシスホールディングス ヤフー リコー 三井住友フィナンシャルグループ 住友大阪セメント 旭硝子 東京海上ホールディングス 野村ホールディングス 電通 => 金融関連

cluster 11: シチズンホールディングス シャープ パナソニック 北越紀州製紙 川崎重工業 日立建機 旭化成

cluster 12: 伊藤忠商事 大阪ガス 日本郵船 日本電気 東京ガス 東急不動産ホールディングス

cluster 13: ANAホールディングス DOWAホールディングス スズキ ファナック 三菱ケミカルホールディングス 住友金属鉱山 東京ドーム

cluster 14: イオン ジーエス・ユアサコーポレーション ソニー ソニーフィナンシャルホールディングス デンカ ニコン ファーストリテイリング 三井不動産 丸井グループ 東ソー 東京建物 東芝 高島屋

cluster 15: 三菱重工業 丸紅 安川電機 日本ハム 日立製作所 本田技研工業 東京エレクトロン 神戸製鋼所

cluster 16: JFEホールディングス アサヒグループホールディングス ダイキン工業 日本板硝子 日本電信電話 日産自動車 東邦亜鉛 豊田通商

cluster 17: T&Dホールディングス あおぞら銀行 セコム 千葉銀行 新生銀行 横浜銀行 第一三共 静岡銀行 => 銀行関連

cluster 18: エーザイ キリンホールディングス フジクラ 三菱倉庫 三菱自動車工業 住友不動産 住友化学 住友商事 住友重機械工業 住友電気工業 協和発酵キリン 大和ハウス工業 大日本印刷 小松製作所 日本たばこ産業 日本電気硝子 明電舎 東京急行電鉄 東洋紡 荏原製作所 鹿島建設 => 住友グループ

cluster 19: NTN パイオニア 三菱電機 京セラ 塩野義製薬 日揮 資生堂

改善点

• 他の情報源からもドキュメントを集める
• 汎用語句の閾値をチューニングする
• トピックの数をチューニングする
• クラスタの数をチューニングする

word2vecで遊ぶ

　MeCabを入れたので、word2vecで日本語文章を解析して遊んでみた。

gensimのインストール
Pythonの場合はgensimをインストールすると、word2vecが使える。
gensimはドキュメントからセマンティックなトピックを自動的に抽出する機能を提供するライブラリ。

pip install --upgrade gensim

使ってみる
太宰治の「人間失格」を使ってモデルを学習させて、単語の足し算、引き算などをやってみた。

# -*- coding: utf-8 -*-

import logging
import codecs
from prettyprint import pp

import MeCab
from gensim.models import word2vec

logging.basicConfig(format='%(asctime)s : %(levelname)s : %(message)s', level=logging.INFO)

with codecs.open('ningen_shikkaku.txt', 'r', 'sjis') as f:
    
    text = f.read()
    tagger = MeCab.Tagger('-Owakati')
    wakati_text = tagger.parse(text.encode('utf-8'))    
    open('wakati.txt', 'w').write(wakati_text)
    
    sentences = word2vec.Text8Corpus("wakati.txt")
    model = word2vec.Word2Vec(sentences, size=100, min_count=1)

    pp(model.most_similar(positive=[u'人生']))
    pp(model.most_similar(positive=[u'人生'], negative=[u'酒']))
    pp(model.most_similar(positive=[u'罪'], negative=[u'罰']))

ちょっとハマったところ

• MeCabはstr型しか受け付けないため、unicode型の文字列を処理したい場合はエンコードしてからMeCabに渡す。
• word2vecでモデルを構築するときにデフォルトだと頻度が5未満の文字は無視されてしまう。頻度がどれくらいなら無視するかはmin_countパラメータで設定できる。テキストにあるはずの ワードをpositive/negativeに入れても「KeyError: u"word xxxxxxx' not in vocabulary"」みたいなエラーが出てハマってしまった。

結果

ドキュメントのサイズが小さかったので、それほど面白い結果はえられなかったが、まあまあ面白かった結果だけ載せておく。

人生 = 哀れ、あさましい、罪、尊い

人生 - 酒 = 受け身、小男、ものしずか

罪 - 罰 = うまい

MeCabをインストール

　遅ればせながらMeCabを使ってみた。

MeCabとは？
形態素解析を行うソフト。
形態素解析とはテキストを最小構成要素に分解し、それぞれの要素の品詞を判別する作業のこと。

インストール
Macにインストールした。

brew install mecab
brew install mecab-ipadic

mecab-ipadicはmecabが形態素解析を行う際に使用する辞書。

使ってみる

$ cat test.txt 
「天は人の上に人を造らず人の下に人を造らず」と言えり。
$ mecab test.txt -o res.txt
「 記号,括弧開,*,*,*,*,「,「,「
天 名詞,一般,*,*,*,*,天,テン,テン
は 助詞,係助詞,*,*,*,*,は,ハ,ワ
人 名詞,一般,*,*,*,*,人,ヒト,ヒト
の 助詞,連体化,*,*,*,*,の,ノ,ノ
上 名詞,非自立,副詞可能,*,*,*,上,ウエ,ウエ
に 助詞,格助詞,一般,*,*,*,に,ニ,ニ
人 名詞,一般,*,*,*,*,人,ヒト,ヒト
を 助詞,格助詞,一般,*,*,*,を,ヲ,ヲ
造ら 動詞,自立,*,*,五段・ラ行,未然形,造る,ツクラ,ツクラ
ず 助動詞,*,*,*,特殊・ヌ,連用ニ接続,ぬ,ズ,ズ
人 名詞,一般,*,*,*,*,人,ヒト,ヒト
の 助詞,連体化,*,*,*,*,の,ノ,ノ
下 名詞,一般,*,*,*,*,下,シタ,シタ
に 助詞,格助詞,一般,*,*,*,に,ニ,ニ
人 名詞,一般,*,*,*,*,人,ヒト,ヒト
を 助詞,格助詞,一般,*,*,*,を,ヲ,ヲ
造ら 動詞,自立,*,*,五段・ラ行,未然形,造る,ツクラ,ツクラ
ず 助動詞,*,*,*,特殊・ヌ,連用ニ接続,ぬ,ズ,ズ
」 記号,括弧閉,*,*,*,*,」,」,」
と 助詞,格助詞,引用,*,*,*,と,ト,ト
言え 動詞,自立,*,*,一段,連用形,言える,イエ,イエ
り 助動詞,*,*,*,文語・リ,基本形,り,リ,リ
。 記号,句点,*,*,*,*,。,。,。
EOS

Pythonから使ってみる
まずはMeCabのPython bindingsをインストール。

pip install https://mecab.googlecode.com/files/mecab-python-0.996.tar.gz

これでPythonからMeCabを使える。
まずは分かち書きをやってみる。

$ python
>>> import MeCab
>>> text = '私はその人を常に先生と呼んでいた。'
>>> tagger = MeCab.Tagger('-Owakati')
>>> print tagger.parse(text)
私 は その 人 を 常に 先生 と 呼ん で い た 。 

名詞だけ取り出してみる。

$ python
>>> import MeCab
>>> 
>>> text = '私はその人を常に先生と呼んでいた。'
>>> tagger = MeCab.Tagger('-Ochasen')
>>> node = tagger.parseToNode(text)
>>> res = []
>>> while node:
...   if node.feature.split(',')[0] == '名詞':
...     res.append(node.surface)
...   node = node.next
... 
>>> print ' '.join(res)
私 人 先生

AWS勉強記（4）機械学習系ライブラリインストール

　EC2（Amazon Linux AMI）にPythonの機械学習系ライブラリを入れた。

• numpy
• scipy
• sklearn

インスタンスタイプはt2.microを選択したが、scipyのビルド中にエラーが出た。

メモリが不足していたようで、swapファイルを作成すると正常にインストールできた。

やったこと

パッケージのアップデート、スワップファイルの作成、ビルドに必要なツールのインストール、機械学習系ライブラリのインストール。

実行コマンドは以下のとおり。

sudo yum update

sudo dd if=/dev/zero of=/swapfile bs=1M count=512
sudo chmod 600 /swapfile
sudo mkswap /swapfile
sudo swapon /swapfile

sudo yum install gcc-c++ python27-devel atlas-sse3-devel lapack-devel

sudo pip install numpy
sudo pip install scipy
sudo pip install scikit-learn

AWS勉強記（3）EC2の起動

EC2とは？

• Elastic Compute Cloudの略
• クラウドでサイズを変更できる仮想コンピューティング環境
• 目的にあわせてインスタンスタイプ（t2.micro、t2.smallなど）を選べる
• インスタンスタイプごとにcpu, memory, storage, gnuの性能/容量が異なる
• 好みのAMI（Amazon Machine Image）を選べる
• AMIはマシンのテンプレート。Linux、windowsなどのOSと追加のソフトウェアをパッケージしたもの。

やったこと

1: http://docs.aws.amazon.com/ja_jp/AWSEC2/latest/UserGuide/get-set-up-for-amazon-ec2.html#create-an-iam-user

• キーペアを作成（sshでログインするときの鍵）
• VPCの作成（AWSリソースが起動する仮想ネットワーク）
• セキュリティグループの作成（ファイアウォールのようなもの）
• port、プロトコルごとに受け付けるIPアドレスを設定できる！

2: http://docs.aws.amazon.com/ja_jp/AWSEC2/latest/UserGuide/ec2-launch-instance_linux.html

• EC2インスタンスの起動

3: http://docs.aws.amazon.com/ja_jp/AWSEC2/latest/UserGuide/ec2-connect-to-instance-linux.html#using-ssh-client

• インスタンスへの接続
• 普通にsshでログインできる
• ~/.ssh/configを書いておくと便利。

ポイント

• インスタンス作成時にインスタンスをどのVPC（仮想ネットワーク）に配置するか決められる
• VPCを決めたら、VPCに定義されているセキュリティグループのうち、どれをインスタンスに適用するか選べる
• インスタンスを停止すると、Public IP、Public DNSが変わってしまう

AWS勉強記（2）IAMでユーザ作成

IAMとは？

• Identity and Access Management
• ユーザの認証と各機能へのアクセス承認を管理するための機能

やったこと

• グループの作成
• ユーザの作成
• ユーザをグループに登録

ポイント

• 一つのAWSアカウントを複数ユーザで共有できる
• AWSを契約したときに作成したユーザはルートユーザ、IAMで作成したユーザをIAMユーザと呼ぶ。
• 通常の用途ではIAMユーザを使う。
• ポリシー（アクセス許可）はグループ単位で設定する。
• IAMユーザ名とパスワードはManagement Consoleにログインするときに使う
• access key id とsecret access keyはAWSのAPIを利用するときに使う

参考リンク

http://docs.aws.amazon.com/ja_jp/IAM/latest/UserGuide/introduction.html

http://docs.aws.amazon.com/ja_jp/IAM/latest/UserGuide/getting-started_create-admin-group.html

http://docs.aws.amazon.com/ja_jp/IAM/latest/UserGuide/getting-started_how-users-sign-in.html

AWS勉強記（1）サインイン

　AWSの勉強を始めたので、やったことを書き溜めていきます。

• 無料お試しでユーザー登録。
• AWSマネジメントコンソールにログイン。
• リージョン（画面の右上から選択）を「東京」に変更。

ウィルコクソンの順位和検定

　2群のデータの差を検定するノンパラメトリック手法。以下のような場面で使う。

• データが順位データの場合（5:とても満足、4:満足、3:普通、2:不満、1:とても不満のようなリッカート尺度など）
• データが量的データだが、正規性を仮定できない場合

例題1

東京と大阪の住人を無作為に8人ずつ選んで、内閣支持に関するアンケートを行った。

5: 支持する、4:やや支持する、3:どちらとも言えない、2:やや支持しない、1:支持しない

以下の結果から、東京と大阪で内閣に対する意見は異なっているか検定せよ。

東京	大阪
4	3
3	3
3	1
2	5
1	2
3	5
4	1

データがリッカート尺度で与えられるため、ウィルコクソンの順位和検定を用いる。

> wilcox.test(c(4,3,3,2,1,3,4), c(3,3,1,5,2,5,1))

 Wilcoxon rank sum test with continuity correction

data:  c(4, 3, 3, 2, 1, 3, 4) and c(3, 3, 1, 5, 2, 5, 1)
W = 25.5, p-value = 0.9475
alternative hypothesis: true location shift is not equal to 0

p値 = 0.9475なので、 帰無仮説は棄却できない。
よって東京・大阪で有意な差があるとは言えない。

例題2
ある新薬の効果を試すために被験者8名を集めた。
まず8名のγ-GTP値を計測し、その後2ヶ月新薬を服用してもらい、2ヶ月後γ-GTP値を再計測した。
結果は以下のとおりである。
新薬に効果があるかどうかを検定せよ。

被験者番号	投薬前	投薬後
1	100	95
2	110	100
3	120	105
4	130	110
5	135	120
6	140	145
7	220	170
8	250	195

まず検定対象の分布に正規性があるかどうか確認する。

> shapiro.test(c(100, 110, 120, 130, 135, 140, 220, 250))

 Shapiro-Wilk normality test

data:  c(100, 110, 120, 130, 135, 140, 220, 250)
W = 0.81118, p-value = 0.03766

すると有意水準5%で投薬前のデータは正規性を持たないことが分かる。

正規性が仮定できないため、パラメトリック手法であるt検定は使うことができない。
よって、ノンパラメトリック手法であるウィルコクソンの順位和検定で検定する。
今回は同じ被験者の投薬前と投薬後の比較をするため、pairedオプションをTrue（群間の対応あり）にする。

> x_before <- c=""> x_after <- c=""> wilcox.test(x_before, x_after, paired=T)

 Wilcoxon signed rank test with continuity correction

data:  x_before and x_after
V = 34.5, p-value = 0.02471
alternative hypothesis: true location shift is not equal to 0

上記の結果から、有意水準5%で新薬に効果ありと言える。

フィッシャーの正確確率検定

　標本サイズが小さい場合（期待度数が5未満のマスが20%以上ある場合）は、ピアソンのカイ二乗検定の検定結果は不正確になる（p値が低く計算されてしまう）ことが知られている。
そのような場合は、フィッシャーの正確確率検定を使う。

例題1
高校生の携帯電話所持状況についてアンケートを行った。
以下のアンケート結果から、男女で携帯の所有状況に差があるかどうかを検定せよ。（有意水準は5%とする）

	持っている	持っていない
男	17	3
女	15	5

> fisher.test(matrix(c(17, 3, 15, 5), ncol=2, byrow=T))

 Fisher's Exact Test for Count Data

data:  matrix(c(17, 3, 15, 5), ncol = 2, byrow = T)
p-value = 0.6948
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  0.3011853 14.0494273
sample estimates:
odds ratio 
  1.859142 

p値が0.6948となるため、 帰無仮説は棄却されない。
よって、男女で有意な差があるとは言えない。

因みに同じ問題をカイ二乗検定で検定すると、以下のように有意差ありと判定されてしまう。Rは以下のように「カイ自乗近似は不正確かもしれません」と警告を出してくれるので親切。

> chisq.test(matrix(c(17, 3, 15, 5), ncol=2, byrow=T), correct=F)

 Pearson's Chi-squared test

data:  matrix(c(17, 3, 15, 5), ncol = 2, byrow = T)
X-squared = 0.625, df = 1, p-value = 0.4292

 警告メッセージ: 
 chisq.test(matrix(c(17, 3, 15, 5), ncol = 2, byrow = T), correct = F) で: 
   カイ自乗近似は不正確かもしれません 

問題2
血液型ごとにお酒が好きか嫌いかのアンケートを行った。
以下のアンケート結果から、血液型によってお酒の好き嫌いに差があるかどうか検定せよ。（有意水準は5%とする）

	好き	嫌い
A型	25	15
O型	23	7
B型	18	2
AB型	5	5

> fisher.test(matrix(c(25,15,23,7,18,2,5,5), ncol=2, byrow=T))

 Fisher's Exact Test for Count Data

data:  matrix(c(25, 15, 23, 7, 18, 2, 5, 5), ncol = 2, byrow = T)
p-value = 0.05008
alternative hypothesis: two.sided

p値 = 0.05008なので、 帰無仮説は棄却されない。
よって血液型に有意な差があるとは言えない。

ピアソンのカイ二乗検定

　ピアソンのカイ二乗検定をRでやってみた。R便利だな。

例題1
あるWebサイトのPV数を曜日ごとにまとめると以下のようになった。
曜日によってPV数に差があるかどうかを検定せよ。（有意水準は5%とする）

曜日	PV数
日曜	400
月曜	591
火曜	669
水曜	704
木曜	731
金曜	528
土曜	382

> chisq.test(c(400, 591, 669, 704, 731, 528, 382))

 Chi-squared test for given probabilities

data:  c(400, 591, 669, 704, 731, 528, 382)
X-squared = 209.9, df = 6, p-value < 2.2e-16

2.2e-16（p値） < 0.05（有意水準）となるため、
「曜日ごとの差はない」という帰無仮説は棄却され、このサイトのPV数は曜日によって差があると言える。

例題2
あるWebサイトにはPC用ページとスマホ用ページが存在する。
それぞれのページにおいて、CV数、非CV数をまとめると下表のようになった。
PC用ページとスマホ用ページでCVに影響があるかどうか検定せよ。（有意水準は5%とする）

	CV数	非CV数
PC用ページ	138	2010
スマホ用ページ	49	952

CV数 / (CV数 + 非CV数)でコンバージョン率を計算すると、
PC用ページ = 0.0642
スマホ用ページ = 0.0490
となり、一見すると有意な差がありそうに見える。

この差が統計的に意味のあるものなのかどうかをピアソンのカイ二乗検定で検定してみる。

> chisq.test(matrix(c(138, 2010, 49, 952), ncol=2, byrow=T))

 Pearson's Chi-squared test with Yates' continuity correction

data:  matrix(c(138, 2010, 49, 952), ncol = 2, byrow = T)
X-squared = 2.5923, df = 1, p-value = 0.1074

0.1074（P値）> 0.05（有意水準）となるため、帰無仮説は棄却されない。
よって、PC用ページ、スマホ用ページで有意な差があるとは言えない。

いろいろな平均

　以下の３つの平均について概要とどのような場面で使うのかをまとめた。
用途に応じて"平均"を使いわける必要があると再認識した。

1. 算術平均（arithmetic mean）
2. 幾何平均（geometric mean）
3. 調和平均（harmonic mean）

算術平均
一般的に使われる平均。観測 値の和を観測値の個数で割った値。

あるバスケットボールチームのメンバーの身長は188cm、187cm、168cm、184cm、197cmであった。メンバーの平均身長を求めよ。

という場合は、(188 + 187 + 168 + 184 + 197) / 5 = 184.8cmが算術平均となる。

幾何平均
観測値の積を(1/観測値の個数) 乗した値。

ある企業の売り上げは
2013年 1.5倍
2014年 2.0倍
2015年 3.0倍
となった。
2013年から2015年の3年間でこの企業の売上は年平均何倍になったか？

という問題を考える。

算術平均を使うと、年平均(1.5 + 2.0 + 3.0) / 3 = 2.167倍の成長となる。
しかしこれだと、3年で2.166666^3 = 10.176倍となってしまい、実際の3年で9倍とは異なる結果となってしまう。

幾何平均を用いると、年平均(1.5 * 2.0 * 3) ^ (1/3) = 2.080倍の成長となる。
これだと、3年で2.080^3 = 9倍となり、実際の3年で9倍と同じになる。

この例のように、成長率、増加率などを扱う場合は、幾何平均を用いるとよい。

調和平均
観測値の逆数の算術平均の逆数。

あなたは、A地点からB地点までを直線で往復した。
往路は時速30キロで起動した。復路は時速50キロで移動した。
往路と復路を平均すると、時速何キロで移動したことになるか？

A地点からB地点までの距離をL kmとする。
往復での移動距離は2Lとなる。往復に要した時間は、L/30 + L/50である。
よって往復での平均速度は
2L / (L/30 + L/50)
= 2/ (1/30 + 1/ 50)
= 1/ ((1/30 + 1/ 50) / 2)
となり、これは、"30の逆数と50の逆数の算術平均"の逆数である。つまり30と50の調和平均である。

この例のように割合（※今回の場合は速度=距離と時間の割合）の平均を扱う場合は、調和平均を用いるとよい。

[TED] Vijay Kumar: The future of flying robots

　ドローンの話。小型ドローンが群れをなして飛びまわる姿はなかなかキモい...

flying robotができること

• ビルの中の様子を撮影しハイレゾな(5cmの解像度)地図を作ってくれる
• 飛行中に物をつかんだり、物を持ったまま飛行したりすることができる
• 小さいロボットほど慣性が小さく、慣性が小さいほど衝突に抵抗がある
• ロボットは隣のロボットを認識することで群れをなすことができる

flying robotの活用例

• 農業への応用
• 果実の個数を数える→収穫量の正確な把握→生産チェーンの最適化
• 葉の相対面積からどれだけ光合成できるか分かる→それぞれの木の健康状態が分かる
• 白化した植物を早期発見
• 収穫量10%増加、原料（水など）25%減少が期待される

二項分布を正規分布で近似できる条件

　二項分布B(n, p)はnが十分に大きい場合、正規分布N(np, np(1-p))で近似できる。
正規分布のパラメータnp, np(1-p)はそれぞれもとの二項分布の平均、分散である。
”nが十分に大きい”というのはアバウトな表現だが、一般にnp>5, n(1-p)>5を満たせば十分に大きいと考えてよい。

　という話を実際にグラフに書いて確かめてみた。

from math import sqrt

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom, norm

p = 0.1
 
# plot binomial distribution
ns = [10, 30, 50]
colors = ['r', 'g', 'b']
for i in range(len(ns)):
    n = ns[i]
    c = colors[i]
    xs = np.arange(0, 20+1)
    ys = binom.pmf(xs, n, p)
    plt.plot(xs, ys, color=c, label="binom n={sample_n}".format(sample_n=n))
 
#plot normal distribution
n = 50
xs = np.linspace(0.0, 20, 100)
ys = norm.pdf(xs, loc=n*p, scale=sqrt(n*p*(1-p)))
plt.plot(xs, ys, color='k', linestyle='--', label="approximate norm n=50")
  
plt.legend()
plt.show()

下が表示されたグラフ。
この場合、np>5, n(1-p)>5を満たすnは50。n=50になると分布の形がほぼ釣鐘型になっているのが分かる。また、黒の点線で表示した正規分布で近似できている。

信頼区間の信頼係数について

　信頼区間の信頼係数（信頼度）が何を意味しているのかよく分かっていなかったので、復習した。

　例えば、信頼係数95%というときの95%が何を意味しているかというと、

1. 母集団からある標本（ex.カブトムシ30匹）を取ってきて信頼区間を求める
2. 母集団から別の標本（ex.上とは別のカブトムシ30匹）を取ってきて信頼区間を求める
3. 母集団からまたまた別の標本（ex.上の2つとは別のカブトムシ30匹）を取ってきて信頼区間を求める

.....（以下同様）

ということを100回行った場合、95回は母数（ex.母集団の平均）が推定した信頼区間に入っている。という意味である。

　これを実際に確認したくて、t分布を用いた母平均の区間推定を行うサンプルプログラムを書いてみた。「平均=172.0、標準偏差=6.0の正規分布に従う母集団から10個の標本をサンプルして、母平均の区間を信頼係数95%で推定する。」を100回行って、母平均が推定区間に含まれる回数が95回程度になることを確認した。

import numpy as np
from scipy import stats
from math import sqrt

def t_interval(xs, alpha):
    n = len(xs)
    t = stats.t.ppf(1.0 - (1.0 - alpha)/2.0, n-1) 
    x_var = np.mean(xs)
    s = np.std(xs)
    return x_var - t * s / sqrt(n-1), x_var + t * s / sqrt(n-1)

mu, sigma = 172.0, 6.0
size = 10

trials = 100
ok = 0
for _ in xrange(trials):
    xs = np.random.normal(mu, sigma, size)
    lb, ub = t_interval(xs, 0.95)
    if lb <= mu and mu <= ub:
        ok = ok + 1

print ok

多変量正規分布の形

　共分散行列の形が変わった時に、多変量正規分布の形がどう変わるのかイメージできなかったのでpyplotで図を書いてみた。

　まず簡単な形から。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

x, y =  np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
pos = np.dstack((x, y))

cov_matrix = [[1,0],[0,1]]
gauss = multivariate_normal(mean=[0,0], cov=cov_matrix)
plt.contourf(x, y, gauss.pdf(pos))
plt.axis('equal')
plt.axis([-3, 3, -3, 3])
plt.show()

これはどういう形になるか予想できる。

次に共分散行列の形を以下のように変えてみる。

cov_matrix = [[1,0],[0,3]]

これも予想できる。y軸方向に長い楕円になる。

次に非対角成分の値を変えてみる。

cov_matrix = [[1,0.5],[0.5,1]]

うーむ、軸がずれるのは予想してたけど、どの向きに軸がずれるのかがイメージできない。ちょっと考えてみて、

_, v = np.linalg.eig(cov_matrix)
print v

とすると、

[[ 0.70710678 -0.70710678]
[ 0.70710678  0.70710678]]

と表示されて、ああそうかと納得した。

固有ベクトルを軸にとってるわけか。
さらに共分散行列は対称行列なので固有ベクトルは直交すると。
よく出来てるなー。

カイ二乗値を計算する

　カイ二乗値を計算して、カイ二乗検定とはそういうことだったのかと今さら納得した。

まず、偏りのないサイコロを振ってカイ二乗値を計算してみる。

import bisect

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2

np.random.seed(1226)

# create a dice with given frequency
def create_dice(freq):
    cdf = np.cumsum(freq)
    def dice():
        p = np.random.uniform()
        return bisect.bisect_left(cdf, p)
    return dice

# calculate chi-squared value
def calc_chi_sq(freq, sample_size):
    dice = create_dice(freq)
    sample = np.bincount([dice() for _ in xrange(sample_size)])
    populatoin = [sample_size * f for f in freq]
    return np.sum([(x-y)**2/y for (x, y) in zip(sample, populatoin)])

freq = [1./6, 1./6, 1./6, 1./6, 1./6, 1./6]

# plot dice histogram
dice = create_dice(freq)
plt.hist([dice() for _ in range(3000)], bins=range(7), normed=True)
plt.show()

# plot chi-squared
chi_sqs = [calc_chi_sq(freq, 3000) for _ in xrange(3000)]
plt.hist(chi_sqs, bins=50, normed=True)
x = np.linspace(0,35,50)
y = chi2.pdf(x, df=5)
plt.plot(x, y, color='r')
plt.show()

実行したら以下の２つのグラフが出力される。
上がサイコロの目の分布。下がカイ二乗値の分布。

カイ二乗値はカイ二乗分布に従っていることが分かる。

続いて、激しく偏った分布を持つサイコロを振って、同じことをしてみる。
上のコードの

freq = [1./6, 1./6, 1./6, 1./6, 1./6, 1./6]

の部分を

freq = [3./12, 1./12, 4./12, 2./12, 1./12, 1./12]

に変えて実行してみると、

なるほど、母集団の分布に関係なくカイ二乗値はカイ二乗分布に従うと。

ということで、カイ二乗検定が何をやってるのか今さらながら納得できた。。。

Spockでテストする

Spcokとは？

• テスティングフレームワーク
• Java、Groovyのアプリケーションのテストができる
• テストはGroovyで書く

導入方法
build.gradleに以下を追記する。

cglib:cglib-nodep:2.2というライブラリを入れておくとインターフェースを作っていないクラスの場合でもモックオブジェクトが作れて便利。

apply plugin: 'java'
+apply plugin: 'groovy'
apply plugin: 'eclipse'

repositories {
    mavenCentral()
}

dependencies {
    compile 'org.projectlombok:lombok:1.16.6'
    compile 'io.netty:netty-all:4.1.0.Beta1'
    testCompile 'junit:junit:4.12'
+    testCompile 'org.spockframework:spock-core:1.0-groovy-2.4'
+    testCompile 'cglib:cglib-nodep:2.2'
}

サンプル1
まずは簡単なやつから。

package com.kenjih.spoc;

public class FizzBuzz {
    public String calc(int x) {
        if (x % 15 == 0)
            return "Fizz Buzz";
        if (x % 3 == 0)
            return "Fizz";
        if (x % 5 == 0)
            return "Buzz";
        return String.valueOf(x);
    }
}

上のFizz Buzzをテストするコードを書いてみる。

package com.kenjih.spoc

import spock.lang.Specification;

class FizzBuzzTest extends Specification {

    def test() {
        setup:
        FizzBuzz f = new FizzBuzz()
        
        expect:
        f.calc(x) == ret
        
        where:
        x      |  ret
        3      |  "Fizz"
        10    |  "Buzz"
        15    |  "Fizz Buzz"
        0      |  "Fizz Buzz"
        101 |  "101"
    }
}

サンプル2
もうちょっと複雑なテストを書いてみる。
ServiceからDaoを呼んでいるときに、Daoの戻り値を恣意的に変えたり、Daoのメソッドが呼ばれた回数をカウントしたりできる。

package com.kenjih.spoc;

import java.util.HashMap;
import java.util.Map;

import lombok.NonNull;
import lombok.RequiredArgsConstructor;

@RequiredArgsConstructor
public class UserService {
    
    @NonNull
    final UserDao userdao;
    
    public boolean login(String uid, String pass) {
        if (uid == null || pass == null)
            return false;
        
        Map<String, Object> cond = new HashMap<String, Object>();
        cond.put("uid", uid);
        cond.put("pass", pass);
        
        return userdao.count(cond) == 1;
    }
}

メソッド名には日本語も使えて便利。。

package com.kenjih.spoc

import spock.lang.Specification;

class UserServiceTest extends Specification {
    def "uidとpassが正しい場合のみログイン成功する"() {
        setup:
        UserDao userDao = Mock()
        userDao.count(_) >> rowCount
        UserService userService = new UserService(userDao)
        
        expect:
        userService.login(uid, pass) == res

        where:
        uid       | pass      | rowCount | res 
        null      | "1234" | 1               | false
        "taro" | null        | 1               | false
        "taro" | "1234" | 0               | false
        "taro" | "1234" | 1               | true
    }
    
    def "uidかpassがnullのときはDBアクセスしない"() {
        setup:
        UserDao userDao = Mock()
        UserService userService = new UserService(userDao)
        
        when:
        userService.login(null, "abcd")
        then:
        0 * userDao.count(_)
        
        when:
        userService.login("jiro", null)
        then:
        0 * userDao.count(_)

        when:
        userService.login("yukiko", "xyz")
        then:
        1 * userDao.count(_)
    }
}