Tech Tips: 12月 2015

2015年12月31日木曜日

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

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

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

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

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

2015年12月28日月曜日

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 パイオニア三菱電機京セラ塩野義製薬日揮資生堂

改善点

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

2015年12月27日日曜日

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"」みたいなエラーが出てハマってしまった。

結果

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

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

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

罪 - 罰 = うまい

2015年12月25日金曜日

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)
私 人 先生

2015年12月23日水曜日

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

2015年12月18日金曜日

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が変わってしまう

2015年12月17日木曜日

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

2015年12月16日水曜日

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

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

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

2015年12月13日日曜日

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

　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用ページ、スマホ用ページで有意な差があるとは言えない。

2015年12月12日土曜日

いろいろな平均

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

算術平均（arithmetic mean）
幾何平均（geometric mean）
調和平均（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の調和平均である。

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

2015年12月9日水曜日

[TED] Vijay Kumar: The future of flying robots

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

flying robotができること

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

flying robotの活用例

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

2015年12月6日日曜日

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

　二項分布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%が何を意味しているかというと、

母集団からある標本（ex.カブトムシ30匹）を取ってきて信頼区間を求める
母集団から別の標本（ex.上とは別のカブトムシ30匹）を取ってきて信頼区間を求める
母集団からまたまた別の標本（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

2015年12月2日水曜日

多変量正規分布の形

　共分散行列の形が変わった時に、多変量正規分布の形がどう変わるのかイメージできなかったので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]]

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

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

2015年12月1日火曜日

カイ二乗値を計算する

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

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

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]

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

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

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

登録: 投稿 (Atom)

Tech Tips

Page List

Search on the blog

2015年12月31日木曜日

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

2015年12月28日月曜日

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

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

2015年12月27日日曜日

word2vecで遊ぶ

2015年12月25日金曜日

MeCabをインストール

2015年12月23日水曜日

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

2015年12月18日金曜日

AWS勉強記（3）EC2の起動

2015年12月17日木曜日

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

2015年12月16日水曜日

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

2015年12月13日日曜日

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

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

ピアソンのカイ二乗検定

2015年12月12日土曜日

いろいろな平均

2015年12月9日水曜日

[TED] Vijay Kumar: The future of flying robots

2015年12月6日日曜日

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

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

2015年12月2日水曜日

多変量正規分布の形

2015年12月1日火曜日

カイ二乗値を計算する

Blogger Syntax Highliter