読者です 読者をやめる 読者になる 読者になる

本日の給油

タイトルに「本日の給油」としてあるけど、ゴールデンウィーク中に行った南東北のカブ旅の給油記録。

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2017-03-18 8559.6 2.76 125.72 45.94 2.74
2017-04-29 8703.9 2.77 127.80 52.09 2.45
2017-04-29 8905.9 3.21 128.97 62.93 2.05
2017-04-30 9032.5 2.20 131.82 57.55 2.29
2017-05-01 9190.3 2.70 127.04 58.44 2.17
2017-05-01 9325.7 2.52 128.97 53.73 2.40
2017-05-01 9469.7 2.17 129.03 66.36 1.94
2017-05-02 9625.3 3.71 133.96 41.94 3.19
2017-05-02 9787.8 2.81 125.98 57.83 2.18

東京から国道6号線を北上して、現時点で原付2種が通れるぎりぎりのところまで進む。高速道路ならば自動2輪車も帰宅困難区域を通れるようになったのだけど110 ccのカブは高速に乗れないので、迂回して仙台に向かった。 f:id:amarui:20170430095027j:plain

下の写真は宮城県の荒浜。後ろに荒浜小学校(Googleマップ)が見える。とにかく鳥の鳴き声がよく聞こえた。 f:id:amarui:20170430153015j:plain

帰りは国道121号を南下。会津鉄道湯野上温泉駅Googleマップ)周辺の桜がきれいだった。 f:id:amarui:20170502091357j:plain

複数箇所のダムを巡ってダムカードをいただいた。ダムではないけど渡良瀬遊水池にもカードがあった。 f:id:amarui:20170502190859j:plain

1100 km超の一般道をガソリン代2857円で走る3泊4日の旅でした。

碁盤の大きさ(19路・13路・9路)

碁盤の大きさとして主なものに何種類かあります。プロの対局にも使われている一般的な大きさは19路盤(じゅうきゅうろばん)ですが、入門用として13路盤と9路盤もよく使われています。他にも6路盤、7路盤、15路盤などもあると聞いたことがありますが、特によく使われているのが9路・13路・19路の三種類なのは、なにか理由があるのでしょうか。

一辺の長さを考えると19路の約半分が9路なので、それらの間の大きさは15路となるのですが、囲碁で使われるのは13路のほうが多いようです。

目の数(≒面積)で考えると、9路は81目、13路は169目、19路は361目と、およそ倍々になっていることが分かります(きちんと計算すると、81→(2.09倍)→169→(2.14倍)→361、となっています)。囲碁は石で囲んだ陣地の広さを競うボードゲームなので、9路盤で入門し、上達にあわせて面積を倍々に増やしていくのが理にかなっているのかもしれません。

19路盤の次の段階があるとすると、361目の倍の722目に近い数として27路盤(27×27=729)でしょうか。9路盤の前の段階として東京大学囲碁講義などの入門で6路盤が使われるのも面積比の関係なのかしら。

Python(とSciPy)で簡単な統計的仮説検定

Excelでもできる程度の簡単な統計的仮説検定ならSciPyの範囲内でできるもよう。以下の実行環境はMac OS X 10.11.6 El Capitan上でHomebrewのpyenvからインストールしたAnaconda 4.3.1(Python 3.6.0、NumPy 1.11.3、Scipy 0.18.1)です。

下記の例は、いずれも伝統的な(あるいは頻度論に基づく)統計的仮説検定をPythonでやる方法をざっくり書いたものですので、統計については別途教科書などを参照して下さい……。

あらかじめ以下のようにscipy.statsモジュールを読み込んでおきます。

>>> from scipy import stats

2項検定

コイン投げのように成功率1/2(=0.5)のもとで、20回中16回表を出すことに成功したとき。

>>> p = stats.binom_test(16, 20, 0.5)
>>> p
0.011817932128906248

有意水準α=.05としたとき、統計的有意となりました。(p値の基準は色々ありますのでそれぞれの分野で適切なαやβやNを用いて判断して下さい。)

対応なしのt検定

24人を12人ずつの2グループに分け、それぞれ試料(例えばケーキの味とかカレーの辛さとか)を5段階評価してもらったとします。その2群の平均値に統計的な有意差があるかを調べます。第1グループのk番目の人と第2グループのk番目の人とが違うので「対応なし」を使います。

>>> d1 = [2,3,1,5,3,4,3,5,4,3,2,3]
>>> d2 = [4,3,4,5,4,4,3,5,4,5,4,4]
>>> t, p = stats.ttest_ind(d1, d2)
>>> t
-2.3213496004719611
>>> p
0.029924340790917867

t値がマイナスの時はd1の平均値<d2の平均値です。効果量とか信頼区間とかも見ると良いと思います。

デフォルトでは分散が等質だとして分析しますが、F検定などで分散が等質でないと出た場合には以下のようにウェルチ検定を行うことができます。

>>> t, p = stats.ttest_ind(d1, d2, equal_var=False)
>>> t
-2.3213496004719611
>>> p
0.032728097618600939

対応ありのt検定

12人がそれぞれ2つの試料の両方を5段階評価してもらった場合には、第1グループのk番目の人と第2グループのk番目の人とが同じになるので「対応あり」を使います。ケーキの評価について考えると、そもそも甘いものが嫌いな人と好きな人などの個人差があります。そういった部分を考慮することができる方法です。

>>> d1 = [2,3,1,5,3,4,3,5,4,3,2,3]
>>> d2 = [4,3,4,5,4,4,3,5,4,5,4,4]
>>> t, p = stats.ttest_rel(d1, d2)
>>> t
-2.9303748521637205
>>> p
0.013681242257293657

1要因の分散分析

2群であればt検定を使いますが、3群以上ではまず分散分析をしてから下位検定をすることが多いようです。1要因の分散分析は以下のように行います。

>>> a = [34, 39, 50, 72, 54, 50, 58, 64, 55, 62]
>>> b = [63, 75, 50, 54, 66, 31, 39, 45, 48, 60]
>>> c = [49, 36, 46, 56, 52, 46, 52, 68, 49, 62]
>>> f, p = stats.f_oneway(a, b, c)
>>> f
0.098615166671485177
>>> p
0.9064161716556407

f比が小さい(p値が大きい)ので、3群の平均値に有意な差はなさそうです。

無相関検定

2変数のあいだの相関が有意であるかを見ます。

>>> d1 = [2,3,4,5,3,4,3,5,4,3,2,3]
>>> d2 = [4,3,4,5,4,4,3,5,6,5,4,4]
>>> r, p = stats.pearsonr(d1, d2)
>>> r
0.50052020128514063
>>> p
0.097450118105491076

rを見ると弱い〜中程度の相関がありそうに思えますが、p値がビミョーです。

その他

scipy.statsモジュールのマニュアルに一覧が載っています。カイ二乗検定、順位和検定、コルモゴロフ−スミルノフ検定、フリードマン検定、などなど盛りだくさん。分布関数や記述統計のプログラムもあるのでいろいろな場面で使えそう。ただ、2要因以上の分散分析が見当たらない。同ページの最後の文は

For many more stat related functions install the software R and the interface package rpy.

(これ以外の統計関係の関数についてはRをインストールしてrpyパッケージを使ってね)

とのことなので、外部プログラムと連携させる必要があるようです。もしRとつなげるのであれば、今はrpyやrpy2よりPypeRのほうが良いというはなしです。ここここに解説があります。

(2要因以上の分散分析はstatsmodelsにも入っているのだけど、まだ繰り返し要因などの実装がないので、現時点ではRを使うしかない感じ。scikit-learnも見てみたい。)

Python再開してみています

15年くらい前には中途半端にPythonを使っていたのですが、最近のデータサイエンス業界でのPythonの盛り上がりが気になったので、Pythonに再入門してみました。(15年前は主としてCとJavaをやっていたのでCythonではなくJythonを使っていましたが、いろいろと不満もあってMatlabに移行したのでした)

ここ1〜2週間ほどPythonSoloLearn.comのウェブサイトとiPhoneアプリで勉強しつつ、これまでMatlabやRで書いたコードをPythonに移植しています。

全てのものが一カ所にまとまっているMatlabと違って、PythonではNumPy、SciPy、Matplotlibなど様々なライブラリに散らばっていて面倒だと感じましたが、とりあえずpyenvでanacondaを入れておけば最小公倍数的に必要そうなものが全部そろいますし普段のMac OS X環境を破壊してしまう心配がないので良いと思います。今後、本腰を入れてPythonを使うかどうかは分かりませんが、そうであれば必要なものだけに絞り込んでインストールしたい気がします。

さて、今日はオーディオファイルを読み込んでスペクトルなどを表示するところまで行きました。以下、撥弦楽器のサンプルを分析したものです。

import soundfile as sf
x, fs = sf.read("sound.wav")

# 自前のツール群を読み込み
exec(open("./audioutils.py").read())

# スペクトルを描画(両線形軸)
plot_spectrum(x, fs)

f:id:amarui:20170323131642p:plain

# スペクトルを描画(両線形軸・範囲指定)
plot_spectrum(x, fs, xlim=(0,5000))

f:id:amarui:20170323132810p:plain

# スペクトルを描画(両対数軸)
plot_magnitude(x, fs)

f:id:amarui:20170323132118p:plain

# スペクトル中心(Hz)の時間変化を描画
sc, t = spectral_centroid_moving(x, fs, window_size=4096, hop_size=1024)
import matplotlib.pyplot as plt
plt.plot(t, sc)
plt.xlabel("Time (sec)")
plt.ylabel("Spectral Centroid (Hz)")
plt.grid()
plt.show()

f:id:amarui:20170323155424p:plain

# 実効値(dB)の時間変化を描画
r, t = rms_moving(x, fs, window_size=4096, hop_size=1024)
plt.plot(t, r)
plt.xlabel("Time (sec)")
plt.ylabel("RMS (dB)")
plt.grid()
plt.show()

f:id:amarui:20170323161423p:plain

plt.show()の代わりにplt.savefig("hoge.pdf")などとすると表示せずファイルに保存。

本日の給油

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2017-03-18 8559.6 2.76 125.72 45.94 2.74

今日は埼玉県南埼玉郡まで。前回給油が去年10月29日に茨城県つくば市に行ったときだったので、それからは近所への買い物くらいしかしてなかったことになる。真夏用のメッシュジャケットしか持っていないというのが原因なのかも? 来冬には冬用のジャケットやグローブ買おうかな。