Mar 24, 2019

Pandasのmelt/pivotとstack/unstackの違い

この間、Pandasを使った教材を流し読みしていると、DataFrameクラスのmeltというメソッドが出てきて、pivotの逆だと説明されていた。それを読んで、meltされたものをpivotして元に戻してみよう思って、すきま時間を集めてのべ1時間くらいがんばったが、できなかった。
ついでに、meltとstackの違いがわからなかった。

今月、1時間くらい連続してがんばったら成功して、meltとstackの違いも理解できたので、ここに控えておく。

以下、コード例はJupyter Notebookで書いたインタラクティブシェル向けの形式のものをそのまま貼り付けており、出力例はJupyter Notebookの出力を加工して作成している。

In [1]:
import pandas as pd
df = pd.DataFrame({'月': ['1月', '2月', '3月'],
                   '京都': [110, 115, 144],
                   '大阪': [263, 283, 309],
                   '奈良': [12, 13, 21],
                  })
df
Out [1]:
京都 大阪 奈良
0 1月 110 263 12
1 2月 115 283 13
2 3月 144 309 21

このようなDataFrameがあるとして、これをmeltして、pivotで元に戻してみる。

In [2]:
df_melted = df.melt(id_vars='月', var_name='地域', value_name='人数')
df_melted
Out [2]:
地域 人数
0 1月 京都 110
1 2月 京都 115
2 3月 京都 144
3 1月 大阪 263
4 2月 大阪 283
5 3月 大阪 309
6 1月 奈良 12
7 2月 奈良 13
8 3月 奈良 21
In [3]:
df_pivoted = df_melted.pivot(index='月', columns='地域', values='人数')
df_pivoted
Out [3]:
地域 京都 大阪 奈良
1月 110 263 12
2月 115 283 13
3月 144 309 21

meltしたものをpivotすると大体戻ったが、「月」がindexになっているのと、meltで加えた「地域」が残っているのが異なるので、修正する。
(ちなみに、reset_indexしただけでは「月」が「地域」に入っておかしなことになる)

In [4]:
df_pivoted = df_pivoted.reset_index()
df_pivoted.columns.name = None
df_pivoted
Out [4]:
京都 大阪 奈良
0 1月 110 263 12
1 2月 115 283 13
2 3月 144 309 21

次に、meltとstackの違いを見てみる。
上のdfをstackすると「月」と地域名が同列に処理されてしまうので、まず「月」をインデックスにしてからstackする。

In [5]:
df2 = df.set_index('月')
df2
Out [5]:
京都 大阪 奈良
1月 110 263 12
2月 115 283 13
3月 144 309 21

この場合はstackすると結果が1列なのでDataFrameでなくSeriesになるので、変数名の先頭をdf_でなくsr_にしている。

In [6]:
sr_stacked = df2.stack()
sr_stacked
Out [6]:
1月 京都 110
大阪 263
奈良 12
2月 京都 115
大阪 283
奈良 13
3月 京都 144
大阪 309
奈良 21

当然ながら、stackしたものをunstackすると元に戻る。

In [7]:
sr_stacked.unstack()
Out [7]:
京都 大阪 奈良
1月 110 263 12
2月 115 283 13
3月 144 309 21

meltでは元々column名だった列(「地域」列)がindexにならなかったが、stackではindexになっている。つまり、meltはcolumn名を新たに加えた通常の列に展開することによって表を変形し、stackはcolumn名をindexの新たな階層に展開することによって表を変形する。このことがmeltとstackの主な違いのようだ。
次のようにして、stackしたもののindexを通常の列に戻すと、meltした結果(Out [2])と同じになる。

In [8]:
df_stacked = sr_stacked.sort_index(level=1).reset_index()
df_stacked.columns = ['月', '地域', '人数']
df_stacked
Out [8]:
地域 人数
0 1月 京都 110
1 2月 京都 115
2 3月 京都 144
3 1月 大阪 263
4 2月 大阪 283
5 3月 大阪 309
6 1月 奈良 12
7 2月 奈良 13
8 3月 奈良 21

See more ...

Posted at 09:38 in PC一般 | WriteBacks (0)
WriteBacks

Mar 09, 2019

ピボットテーブルを理解した

15年くらい前にExcelのピボットテーブルを使ってみて以来ずっと、ピボットテーブルとは何なのかが理解できなかった。2年くらい前にも、ピボットテーブルを初心者向けに詳しく解説したITProの記事を読んで、今度こそ理解しようと意気込んで、何も見ずに演習問題を解けるようにもなったが、結局「ピボットテーブル」の意味を理解できず、ただ演習問題の解き方を丸暗記したような感じで終わってしまった。
筆者にとってピボットテーブルは、わかりそうなのにわからない、何故分からないのかがわからない、まるで眼の盲点のように、筆者の脳はそれが理解できない作りになっているかのように思えるものだった。

最近、久々にPandasを勉強する機会があって、ピボットテーブルが出てきて、また理解できなかった。しかし、何回かgroupbyの練習をした後、ふと、ピボットテーブルってgroupbyと似てるなと思い、Webで調べたらピボットテーブルはGroupByを2次元にしたようなものというような説明が見つかって、遂に理解に成功した。

GroupByの動作は、split-apply-combineと言われるように、データをキーによって分割し、分割した単位で何らかの処理を適用し、結合して1つのデータにするものである。
次の図は、データをkey列の値によって分割し、分割された組毎に合計値を計算して1つの値に「集約(aggregate)」し、1つのデータにまとめる、というGroupByの処理の例を表している。

●コード例
import pandas as pd
df = pd.DataFrame({
    'key': ['A', 'B', 'C'] * 3,
    'val': range(1, 10)})
print(df)
print(df.groupby('key').aggregate(sum))
●実行結果
  key  val
0   A    1
1   B    2
2   C    3
3   A    4
4   B    5
5   C    6
6   A    7
7   B    8
8   C    9
     val
key     
A     12
B     15
C     18

Applyは集約に限らず、単に表を変形したり、グループ毎の値を使って値を変換することなども含まれるが、グループ毎に集約するのがGroupByやピボットテーブルの醍醐味だと思うので、この記事の例では集約にしている。

ピボットテーブルは、split-apply-combineを行方向と列方向に同時に行うものと捉えることができる。
次の図は、データをkey1列の値によって行方向に、key2列の値によって列方向に分割し、同様に分割毎に合計値に集約し、1つのデータに結合する、というピボットテーブルの処理の例を表している。

●コード例(pivot_tableメソッド使用)
import pandas as pd
df = pd.DataFrame({
    'key1': ['A', 'B'] * 6,
    'key2': ['①', '②', '③'] * 4,
    'val': range(1, 13)})
print(df)
print(df.pivot_table(index='key1', columns='key2', values='val', aggfunc='sum'))
●実行結果
   key1 key2  val
0     A    ①    1
1     B    ②    2
2     A    ③    3
3     B    ①    4
4     A    ②    5
5     B    ③    6
6     A    ①    7
7     B    ②    8
8     A    ③    9
9     B    ①   10
10    A    ②   11
11    B    ③   12
key2   ①   ②   ③
key1            
A      8  16  12
B     14  10  18

ピボットテーブルはGroupByを2次元にしたようなもの、ということを確認する為に、pivot_tableの代わりにgroupbyを2つ使って同じ結果にしてみる。

●コード例(groupbyメソッドを2つ使用)
print(df.groupby('key1').apply(lambda x: x.groupby('key2')['val'].apply(sum)))
●実行結果
key2   ①   ②   ③
key1            
A      8  16  12
B     14  10  18

See more ...

Posted at 20:26 in PC一般 | WriteBacks (0)
WriteBacks

Dec 13, 2018

Raspberry Pi 2にscikit-learn 0.20をインストール

aptでpython-sklearnをインストールして、pythonで

from sklearn.neural_network import MLPClassifier

すると、見つからないというエラーになった。
Raspbian 8.0 (jessie)の現在最新のaptのpython-sklearnのバージョンは0.14.1-3だが、sklearn.neural_network.MLPClassifierはscikit-learn 0.18以降にしか無いようだ。

そこで、pipで新しいバージョンのscikit-learnをインストールしようとして、

sudo pip install --upgrade pip

するとpip3が壊れたりしたので、
sudo pip uninstall pip
sudo apt install --reinstall python-pip python3-pip

として修復し、virtualenvを使うことにした。

aptでpython-virtualenvをインストールして、

virtualenv (workdir)
cd (workdir)
. bin/activate
pip install --upgrade pip
pip install scikit-learn

とすると、scipyのmake中にblasやlapackが無いというエラーになった。
https://scikit-learn.org/stable/developers/advanced_installation.html
によると、scikit-learnのインストールにはlibatlas-devとlibatlas3-baseが必要とのことだが、scipyのmakeには他に、liblapack-devか、もしくはatlasに合わせるならlibatlas-base-devが必要のようだ(どちらでも成功した)。

aptでlibatlas-dev, libatlas3-base, libatlas-base-devをインストールして、改めて

pip install scikit-learn

すると、数十分後にRaspberry Pi 2からの応答が無くなった。scipyのmake中にメモリのスワップ領域が枯渇するようだった。

スワップ領域を1GB追加しても不足したので、2GB追加した。

su -
dd if=/dev/zero of=/swapfile bs=1024 count=2097152
chmod 0600 /swapfile
mkswap /swapfile
swapon /swapfile

これによって、数時間かかるが、scikit-learnのインストールができるようだ。

ただ、時間がかかるし、makeが失敗して時間を取られるリスクがあるので、やっぱりaptでインストールできるものはaptでインストールすることにした。MLPClassifierと合わせてmatplotlibを使いたいが、もはやpipでmatplotlibをインストールする気にはならない。
現時点で手っ取り早くscikit-learn 0.20.xをインストールするには、以下のようにすれば良いと思う。

sudo apt install python-numpy python-scipy python-virtualenv libatlas-dev libatlas3-base libatlas-base-dev
virtualenv --system-site-packages (dir)
cd (dir)
. bin/activate
pip install --upgrade pip
pip install scikit-learn

※(dir)は任意のディレクトリ名

See more ...

Posted at 23:29 in PC一般 | WriteBacks (0)
WriteBacks

Dec 10, 2018

Raspberry Pi 2にUSBカメラを接続

昨年辺りから、Raspberry Piにカメラを付けて、Deep Learningさせたニューラルネットワークを使って物体認識させるのが流行ってるというか、今年の夏にTensorFlowが公式にRaspberry Piをサポートしたことで、定番から基本になりつつあるようだ。
近い内にAIを勉強しようと思っており(思い始めてからだいぶ経つが)、かつ、Raspberry Piを所有している身としては、一度やっておかないといけないと思っていた。

しかし、Raspberry Piのカメラモジュールは4,000円くらいするので、二の足を踏んでいた。筆者は貧乏性で、確実に元が取れると思えない限りはお金を使わない、よくそれで損をするタイプなのである。
そんな折、昔Skype英会話用に買ったUSBカメラが10回くらい使っただけで眠っているのを思い出した。2010年に買った、Logicool HD Webcam C270というやつである。思い出してからも、どうせ消費電力的に無理だろうと思っていたが、先月ふと思い立って繋いだら動いて、今の所特に電源異常は起こっていない。少し調べた限りでは、このカメラの消費電力が1.0WくらいでRaspberry Pi 2の消費電力が最大4.5W、使用しているACアダプターの出力が9.0Wなので、大丈夫のようだ。

PythonとCでOpenCVを使って動かしたが、安定的に動かすまでには結構試行錯誤したので、それなりに動くようになったコードを貼っておく。

●Pythonのサンプルコード

#!/usr/bin/python
import cv2

WINNAME = "Capture"
FRAME_INTERVAL = 30 # msec

CAPTURE_WIDTH = 1280
CAPTURE_HEIGHT = 720

ret = False
while ret == False:
	cap = cv2.VideoCapture(0)
	cap.set(cv2.cv.CV_CAP_PROP_FRAME_WIDTH, CAPTURE_WIDTH)
	cap.set(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT, CAPTURE_HEIGHT)
	ret, img = cap.read()

key = 0
while key != ord('q'):
	ret, img = cap.read()
	cv2.imshow(WINNAME, img)
	cv2.moveWindow(WINNAME, 0, 0)
	key = cv2.waitKey(FRAME_INTERVAL)
	if key == ord('s'):
		cv2.imwrite('capture.jpg', img)

cap.release()
cv2.destroyAllWindows()

実行する前に、

sudo apt install python-opencv
などとしてaptでpython-opencvをインストールしておく。

なお、Python3でも動かそうと思ったが、aptにpython3-opencvが無く、代わりにaptでpython3-pipをインストールして

sudo pip3 install opencv-python
とすれば良いと思うのだが、下のようなエラーになったので、今回は諦めた。
(pip3 install --upgrade pipとしてpip3を最新にしても同じ結果だった)
Could not find any downloads that satisfy the requirement opencv-python
Cleaning up...
No distributions at all found for opencv-python

●Cのサンプルコード

#include <cv.h>
#include <highgui.h>

#define WINNAME "Capture"
#define FRAME_INTERVAL 30 /* msec. */

#define CAPTURE_WIDTH 1280
#define CAPTURE_HEIGHT 720

int main(int argc, char *argv[])
{
	CvCapture *capture = NULL;
	IplImage *frame = NULL;
	int key;

	while (!capture)
	    capture = cvCreateCameraCapture(0);

	cvSetCaptureProperty(capture, CV_CAP_PROP_FRAME_WIDTH, CAPTURE_WIDTH);
	cvSetCaptureProperty(capture, CV_CAP_PROP_FRAME_HEIGHT, CAPTURE_HEIGHT);

	cvNamedWindow(WINNAME, CV_WINDOW_AUTOSIZE);
	cvMoveWindow(WINNAME, 0, 0);

	do {
		frame = cvQueryFrame(capture);
		cvShowImage(WINNAME, frame);
		key = cvWaitKey(FRAME_INTERVAL);
		if (key == 's')
			cvSaveImage("capture.jpg", frame, NULL);
	} while (key != 'q');

	cvDestroyWindow (WINNAME);
	cvReleaseCapture(&capture);
	return 0;
}

コンパイルする前に、

sudo apt install libopencv-dev
などとしてaptでlibopencv-devをインストールしておく。

コンパイルは、次のコマンドで行った。

gcc camera_test.c `pkg-config --cflags --libs opencv` -lm -o camera_test

いずれも、sを押すとカメラ画像がcapture.jpgに保存され、qを押すと終了する。

See more ...

Posted at 17:36 in PC一般 | WriteBacks (0)
WriteBacks

Nov 28, 2018

MovableTypeのログイン状態が続かなくなった

ある時から急に、このWeblogを管理するMovableTypeの画面で何かする度に、"Invalid Login."と表示されたログイン画面に戻るようになってしまった。WebブラウザのCookieが無効の時のような動作である。
それでもブラウザにUsernameとPasswordを自動入力させながら少しずつ進むならまだ良かったのだが、悪いことにエントリー(記事)の保存ができなかった。エントリーの入力画面で"Save"ボタンを押すと"Invalid Login."と書かれたログイン画面になり、UsernameとPasswordを入力してログインすると、エントリーが保存された時の画面でなく、ログインして最初に開く画面に戻るのだった。

別のブラウザでも同様、Apache2を再起動しても同様だった。
最初にApache2のCookieの動作を疑ったせいか、Webで検索しても同じような症状が書かれたページがなかなか見つけられず、何時間も失ってしまった。
結果的に、MySQLデータベースの'mt_session'テーブルが破損したのが原因だとわかり、'mt_session'で検索すると、'mt_session'テーブルが破損してMovableTypeがこの状態になったという情報がたくさん見つかった。

mysqlコマンドでMySQL DBに接続して、describe mt_session;とすると、

Table 'mt_session' is marked as crashed and last (automatic?) repair failed

と表示され、バックアップをリストアしないといけないくらいに大変なことになったかと思ったが、幸いなことに、
repair table mt_session;

とすると一発で直った。

ついでに、mt_sessionに大量のUser Session(US)の行が残っていたので、UNIXのシェルで

remove_old_sessions --ttl 30 --kind US

を実行して削除した。

See more ...

WriteBacks

Jul 15, 2018

大局観とコンピューターの評価とのずれ

昨日は毎年参加している、地域の将棋大会だった。
今年は家庭が多忙で準備不足だったので、良くて1回戦敗退だろうなと期待せずに参加したら、予想通り、危なっかしく予選突破して、決勝トーナメントの1回戦で敗退した。

予選は2勝1敗、2勝は同じ人から(4人のグループで2戦目は勝った者同士、負けた者同士で対局するので1/3の確率で3戦目は初戦と同じ組み合わせになる)で、2局とも角換わり棒銀で一時必敗になりながら逆転勝ちした。筆者はネット将棋では角換わり棒銀で勝てることが少ないが、この大会の予選ではほぼ全勝で相性が良いようである。

決勝の1回戦の相手は昨年に続きまた優勝常連のI藤さんだった。その時点で準備万端だろうが体調万全だろうが関係なく1回戦敗退が決定したのだが、昨年に続き、筆者に感触の良い手が出た。


この△5七銀である。(筆者が後手)

しかし、昨年と違って、後で激指14で解析すると、この局面の評価値はほぼ±0で、それほど良くなかった。
筆者はこの△5七銀でかなり有利になったと思っていたのだが、そうでもないのか。

しかも、予定通りに飛車を成り込んだ、

この局面では少しマイナスの値(後手有利)に傾いたものの、ここから▲4七銀打 △3七角成 ▲同桂と進んだ局面はプラスの値(先手有利)になっている。

上記の△5七銀の局面を目指して、結構頭を使ってその局面に持ち込んだのだが、無駄な努力だったようだ。このように大局観というか局面の評価が間違ってるのが、筆者がアマ三段から上に行けない原因なのだろう。
確かに、全て清算した後は後手の駒損だし、後手が飛車を成り込んでも先手の方が陣形が良く、持ち駒の銀で固めることができて、容易に攻めれそうにない。
それでも、筆者としては、飛車を成り込むのはとても気持ち良かったのである。

その後、相手がI藤さんであることを敢えて意識しないようにして、気持ち良く指し続けて、次の局面になった。

この局面で▲5三角打とされ、△3四竜と桂馬を抜いて、必勝になったと思ったら、次の▲4四歩でしびれてしまった(と思った)。

局後の感想戦でも、△3四竜は見落としで、一瞬投げようかと思った、とI藤さんが言ってたし、横で見ていた知人と予選で2局指した人も、△3四竜で筆者の勝ちになったと思ったが次の▲4四歩でわからなくなったと言っていた。
しかし、激指14の評価は△3四竜の時点で+500(先手有利)、▲4四歩でも特に変わらず(しかも▲4四歩は最善ではない、最善は▲3一角成〜▲5三角成)、△5二桂▲4三歩成△同竜の次の▲4四金が疑問手で評価値が-200に振れ、次の△4一竜が悪手で+1000に振れていた。

筆者は▲4四歩の時点で必ずこの▲4四金まで進むと思っていたし、実際この▲4四金でどうしようも無くなったと思ったが、これが疑問手とは驚いた。ここで△5三竜 ▲同金 △5八歩成とすれば互角だったようだ。

一局を通して筆者が優勢だった区間はほとんど無かったようだ。どこで間違ったのだろう、と思って解析したら最初から全て間違っていたということであり、筆者の大局観はそれほどまでに間違っているのかと驚かされた。

昨年はこれほどまでには激指14の評価と自分の評価がずれていなかったと思う。筆者はここ半年くらい、パソコンで将棋ウォーズの3分将棋をよくやっているが、3分将棋は苦手でずっと初段である(10分将棋と10秒将棋は三段)。すると対戦相手も初段レベル前後なので、大局観がそのレベルで勝ちやすいかどうかになってしまったのかも知れない。

See more ...

Posted at 15:54 in 将棋 | WriteBacks (0)
WriteBacks

Apr 30, 2018

火曜日生まれの男の子が...

子供を2人持つ母親に、「火曜日生まれの男の子がいますか?」と聞いたら、「はい」と答えた。もう1人の子供も男の子である確率は?

2010年7月、筆者が数学好きだということを知っていた職場の後輩から、帰り道にこんな問題を、答えは1/2でないという注釈付きで教わった。それに対して、筆者が夜中に結構悩んでやっと問題の意味を理解し、13/27だと計算したと送ったメールを先週発掘したので、今は亡き後輩を偲んでもう一度夜中にやってみたら、結構悩んだ。

例えば、上記の「火曜日生まれ」を取り除いて、問題を

子供を2人持つ母親に、「男の子がいますか?」と聞いたら、「はい」と答えた。もう1人の子供も男の子である確率は?
に変えると、問題の意味がわかってくる。2人の子供のどちらか一方について「男の子ですか?」と聞いたのではないので、もしどちらか一方だけが男の子であれば、この「はい」によって既に答えられているのである。「男の子がいますか?」に「はい」と答えられる確率は1/2ではなく3/4なので、その分、その中でもう1人も男の子である確率は1/2より低くなる。
従って、問題の意味は「少なくとも1人が男の子である」場合の中の「2人とも男の子である」確率であり、答えは(2人とも男の子である確率)÷(少なくとも1人が男の子である確率)= (1/4)÷(3/4) = 1/3ということになる。
図で書くと、次のように表せるだろうか。
少なくとも1人男の子がいるのは黄色部分3パターンであり、その中でもう一人も男の子なのは*部分1パターンだけなので、その確率は1/3である。
つまり、この問題も冒頭の問題も、母集団を限定した確率を問う問題、即ち条件付き確率の問題である。

従って、冒頭の問題の答えは、(1人が火曜生まれの男の子、もう1人も男の子の確率)÷(少なくとも1人が火曜生まれの男の子である確率) = ((1/14)*(1/2) + (1/2)*(1/14) - (1/14)2) / (1 - (1 - 1/14)2) = 13/27 となる。
同じように図で書くと、次のようになる。















男日
男月
男火
男水
男木
男金
男土
女日
女月
女火
女水
女木
女金
女土
母集団が黄色の部分に限定された場合の中の*の割合なので、13/27である。

すっきりして、きっと有名な問題なんだろうなと思って調べてみたら、ネタ元はBBC Newsの2010年6月の記事、"Tuesday boy"だとわかり、次のような問題文が書かれていた。

"I have two children. One is a boy born on a Tuesday. What is the probability I have two boys?"
記事の初めにこれの答えは13/27と書いてあるが、見た瞬間、この問題文だと答えは13/27にはならないのではないかと思った。
その違和感が何に起因するのかがはっきりとわからず、1週間悩んでいる。

冒頭の問題で、「はい」と答えた母親ばかり集めれば、確率は13/27になることは、例えば次のようなプログラム(in Python)でシミュレーションすれは確認できる。

from random import randint

N = 100000
n = 0
hit = 0
while n < N:
    boy1 = randint(0, 1)	# 1: boy, 0: girl
    day1 = randint(1, 7)	# 1: Sun, 2: Mon, 3: Tue, ...
    boy2 = randint(0, 1)
    day2 = randint(1, 7)
    if boy1 == 1 and day1 == 3:
        n += 1
        if boy2 == 1:
            hit += 1
    elif boy2 == 1 and day2 == 3:
        n += 1
        if boy1 == 1:
            hit += 1

print(hit / N)
ではBBC Newsの問題はこのシミュレーションに当てはまるだろうか?と考えると、筆者には当てはまるように思える時もあるし、そうでないように思える時もある。

筆者にとって、BBC Newsの問題の違和感の原因の1つは、"born on a Tuesday"といった情報が親から自発的に付け加えられていることである。これによってもう一方の子が男の子である確率が1/3から変化するのであれば、親がその"a boy"に関する情報をさらに付け加えればまたもう一方の子に関する確率が変化することになる。

"One is a boy"に情報を付け加えれば付け加えるほど、重複する確率が下がり、求める確率は1/2に近づく。「少なくとも一方が」の意味が薄れ、もう一方の子の情報が無くなっていくとも言えよう。そうであれば、用いるべきは最も求める確率を左右する情報である"One is a boy"のみではないか、従って求める確率は1/3ではないかと思えてくる。
逆に、"One is a boy"も自発的に出された情報なので、「少なくとも1人は」男の子という意味で出されたのでない、例えば年上の子など、1人だけについて言ったものかも知れないと考えれば、もう1人については何も情報が無く、答えは1/2のように思える。

それに対し、冒頭の問題では水曜日生まれの男の子がいても「はい」とは答えられなかったので、「火曜日生まれ」という情報を無視することはできないし、どちらの子かを特定せずに質問したので、「少なくとも1人は」火曜日生まれの男の子だということになり、条件付き確率は間違いなく13/27である。

See more ...

Posted at 23:15 in 数学 | WriteBacks (0)
WriteBacks

Apr 08, 2018

データサイズが異なる処理時間の差の検定を考える

あるソフトウェアをバージョンアップしたら、処理時間が少し伸びた。
計測結果は、次のような感じだった。

入力データ処理時間(before)処理時間(after)
1300360
2200180
3100120
4500480
5400440
612001300
7250270
8200210
915001800
1010001100

ぱっと見で明らかに悪化しているが、念の為、統計学的に有意な差であることを示しておこうと思って、Excelの分析ツールの「t検定: 一対の標本による平均の検定」(いわゆる対応のあるt検定)を使ったら、p値が意外と小さくなかった。例えば上記のデータだと両側検定で0.0706となる。
それでも、有意水準を10%とすると有意差ありと言えないことはないので、不本意ながらその結果を付けて報告した。

そのメールを発信した日の帰り道に、このt検定は妥当な計算がなされていないのではないかと疑問に思った。 同じ入力データのセットを使ったのでbeforeとafterで対応があるが、入力データはサイズがまちまちなので、(after-before)の差それぞれは同列に扱えない。(after-before)の単純な差でなく、差分の比率、それぞれ何%増であるかが検定の対象である。 専門的に言うと、t検定は標本の母集団が正規分布に従うのが前提であるが、この(after-before)の差は正規分布に従うはずがないので、t検定を用いるのは不適切だということになるだろう。

もしかして、Excelの分析ツールの「t検定: 一対の標本による平均の検定」はそれも考慮したp値を計算してたりしないだろうか、と期待したが、結果を確かめてみると、一般的な「対応のあるt検定」と同様、
t=\frac{\bar{d}}{s/\sqrt{n}}
(dは対応のある数値の差、dの上のバーは平均、sは不偏標準偏差、nは標本データ数)
というt値と、それに対応する自由度n-1のt分布における外側確率だった。

こういった、差分の比率を検定したいケースはよくあると思うので、何か定番の方法があるのではないかと思って少し探したが、全く見つからなかった。

対応するデータの差分を、正規分布に従うように補正(正規化)してt検定すれば良いのだが、よく考えると、一般的には差分をbeforeの値で割ったものが正規分布に従う保証は無いのである。今回のデータが、計算量が入力データのサイズに比例する処理時間であることがわかっているから、差分をbeforeで割ったら大よそ正規分布に従うという仮定を置けるが、もし計算量のオーダーがO(n)でなくO(n2)だったら、差分のオーダーもO(n2)ならbeforeの平方根で割るべきかも知れないし、差分のオーダーが不明なら正規化のしようが無いかも知れない。

今回は結局、対応するデータの差をbeforeとafterの平均で割ることによって正規化し、t検定をやり直した。単純にbeforeの値で割って正規化しても良いだろうが、before→afterとafter→beforeを対称に扱う為、例えば80→100と100→80が相殺されるようにする為にこのようにした。考え過ぎか。
例えば上記のデータだと、このように正規化して計算すると、両側検定のp値が0.0267となる。

See more ...

Posted at 16:31 in 数学 | WriteBacks (0)
WriteBacks

Mar 19, 2018

ディリクレ分布を直観的に理解する

サイコロを振って1,2,3,4,5,6の目が出た回数がそれぞれn1,n2,n3,n4,n5,n6だった時、それぞれの目が出る確率p1,p2,p3,p4,p5,p6は(n1,n2,n3,n4,n5,n6)をパラメーターとするディレクトリ分布に従う、という話を聞いたので、ちょっと勉強しておくことにした。

コインの表が出る確率がpの時、N回振って表が出る回数nはNとpをパラメーターとする二項分布に従う。
コインをN回振って表が出た回数がn回だと観測された時、点推定では表が出る確率はn/Nとなるが、ベイズ推定では表が出る確率(事後確率)pは(n+1, N-n+1)をパラメーターとするベータ分布に従う。
勝手な表記だが、回数がnになる確率P(n)が
P(n;N,p)=Binomial\_PMF(n;N,p)=\pmatrix{N \cr n}p^n(1-p)^{N-n}
(PMF=Probability mass function)という感じであり、事後確率がpになる確率の密度関数f(p)が
f(p;N,n)=Beta\_PDF(p;n+1,N-n+1)=\frac{\Gamma(N+1)}{\Gamma(n+1)\Gamma(N-n+1)}p^n(1-p)^{N-n}=\pmatrix{N \cr n}p^n(1-p)^{N-n}
(PDF: probability density function)ということになる。
この関係のことを、ベイズ推定の用語で、尤度関数が二項分布の場合、ベータ分布は共役事前分布(conjugate prior)であると言うらしい。

多項分布とディリクレ分布もその関係にある。
サイコロをN回振って1〜6の目が出る確率がp1〜p6の時、1〜6の目が出る回数n1〜n6はp1〜p6をパラメーターとする多項分布に従う。(n1+n2+n3+n4+n5+n6=N, p1+p2+p3+p4+p5+p6=1)
サイコロをN回振って1〜6の目が出た回数がそれぞれn1〜n6回だと観測された時、推定される1〜6の目が出る事後確率p1〜p6は(n1+1, n2+1, n3+1, n4+1, n5+1, n6+1)をパラメーターとするディリクレ分布に従う。
nをn1〜n6のベクトル、pをp1〜p6のベクトルとすると、回数がnになる確率P(n)が
P(n;p)=Multinomial\_PMF(n;p)=\frac{\left(\sum_{k=1}^{6}n_k\right)!}{\prod_{k=1}^{6}n_k!}\prod_{k=1}^{6}p_k^{n_k}
であり、事後確率がpになる確率の密度関数f(p)が
f(p;n)=Dirichlet\_PDF(p;n+1)=\frac{\Gamma\left((\sum_{k=1}^6n_k)+1\right)}{\prod_{k=1}^6\Gamma(n_k+1)}\prod_{k=1}^6p_k^{n_k}=\frac{\left(\sum_{k=1}^6n_k\right)!}{\prod_{k=1}^6n_k!}\prod_{k=1}^6p_k^{n_k}
ということである。

ディリクレ分布はベータ分布の確率変数の次元を拡張したものである。
上記のコインの例も多項分布とディリクレ分布で表現でき、コインの表が出る確率がp1、裏が出る確率がp2の時、N回振って表が出る回数n1と裏が出る確率n2は(p1, p2)をパラメーターとする多項分布に従う。
コインをN回振って表が出た回数n1回、裏が出た回数がn2だと観測された時、表が出る事後確率p1、裏が出る事後確率p2は(n1+1, n2+1)をパラメーターとするディリクレ分布に従う。

図1はパラメーターを(3n+1, 2n+1)、n=1〜10としたベータ分布のグラフである。 コインを投げて5回中3回、10回中6回、15回中9回、...、50回中30回、表が出た時の表が出る確率の分布に対応する。
図1
5回中3回でも表が出る確率が0.6である確率が最も高いが、他の確率である確率もそれなりに高いのに対し、回数を増すほど確率が0.6近辺に限定されていく様子がわかる。
最尤推定のような点推定では5回中3回でも50回中30回でも単に0.6であり、その尤もらしさが区別されない。

図2は同様にパラメーターを(3n+1, n+1)、n=1〜10としたベータ分布のグラフである。 コインを投げて4回中3回、8回中6回、...、40回中30回、表が出た時の表が出る確率の分布に対応する。 図2

図3はパラメーターを(2+1, 3+1, 5+1)とした3次元のディリクレ分布のグラフである。 3面しか無いサイコロを10回振って各面が2回、3回、5回出た時の各面が出る確率の分布に相当する。
図3
ちょっとややこしいが、三角形の各頂点がサイコロの各面に対応し、三角形内の各頂点への近さが各面の出る確率に対応し、上方向の高さがその確率の組み合わせになる確率であり、高いほど高温の色になっている。右下の頂点が1つ目の面、右上の頂点が2つ目の面、左の頂点が3つ目の面に対応する。XYZ空間で右下の頂点が(1,0,0)、右上の頂点が(0,1,0)、左の頂点が(0,0,1)にあるとすれば、三角形内の点のX座標、Y座標、Z座標が各面の出る確率である。
上から見ると、図4のようなヒートマップになる。
図4
(0.2, 0.3, 0.5)に対応しそうな所が頂点になっている。

図5はパラメーターを(20+1, 30+1, 50+1)とした3次元のディリクレ分布のグラフである。 3面しか無いサイコロを100回振って各面が20回、30回、50回出た時の各面が出る確率の分布に相当する。
図5
図6は上から見た図である。
図6
それぞれの面が出た割合が同じでも、回数を重ねた方が確率が高い範囲が限定される様子がわかる。

See more ...

Posted at 18:17 in 数学 | WriteBacks (0)
WriteBacks

Feb 01, 2018

(Emacs)Jediで数字に不正な補完候補が出る

筆者はEmacsでPythonのコードを書く時に補完機能としてJedi+CompanyModeを使っているが、筆者の環境では数字に対しておかしな補完が働く。

例えば、

このように、数字を打つと"and"などが補完候補として現れ、Enterキーを押すと

このように、打った数字が消えて"and"になってしまう。

そして慌てて"and"を消してもう一度数字を入力してEnterを押すと、また"and"になるのである。
これは不便である。

CompanyModeでなくAutoCompleteだと同じことにはならないが、やはりおかしな補完が起こることがある。

このように、数字を1文字打つと"if"などが補完候補として現れ、Enterキーを押すと

このように、"if"が足される。

筆者の環境は以下である。
・macOS Sierra バージョン10.12.6
・Emacs 24.5.1 (https://emacsformacosx.comからダウンロード)
 ※同サイトからダウンロードしたEmacs 25.2でも起こる
・jedi 20160425.2156
・company-jedi 20151216.1921

.emacsは次のものだけにしても再現する。

(require 'package)
(package-initialize)

;; for Jedi + company
(require 'company)
(add-hook 'python-mode-hook 'company-mode)
(add-to-list 'company-backends 'company-jedi)

;; for Jedi + auto-complete / company-mode
(require 'jedi-core)
(setq jedi:complete-on-dot t)
(add-hook 'python-mode-hook 'jedi:setup)

Webで調べても、同じ症状の話はほとんど見つからない。
唯一見つけたのは、https://github.com/jorgenschaefer/elpy/issues/1115である。昨年に見つけて時々チェックしているが、今年に入って、筆者と同じ症状の報告が追加された。しかし、私には再現できていないという書き込みや、Jediの問題ではない、Python的に間違った補完ではない(??)、Elpyで回避可能といった書き込みがあり、すぐに修正されそうには見えない。

とりあえずの対策として .emacs に

(setq company-idle-delay 0.5)

を加えて、キー入力から0.5秒間は補完が働かないようにした。

See more ...

Posted at 23:10 in PC一般 | WriteBacks (0)
WriteBacks