統計・確率のお勉強

統計学を中心に色々勉強するブログ

<時系列データ解析> 株価をプロット〜Pythonと株価データを使って勉強する〜

主に参考にする本

「現場ですぐ使える時系列データ分析 データサイエンティストのための基礎知識」を中心に勉強して行きたいと思う。
この本にした理由は

  • 内容がそこまで難しくない
  • Rで実行し、手を動かしながら学べるスタイル
  • 扱うデータが株価と自分の興味と合致している

の3つある。私は普段はRをほとんど使わないで言語はPython(Anaconda)、環境はJupyter notebookを使っているので、この本に載っている内容をPythonでやって行きたいと思う。

株価をプロットする

何を分析するにしても、まずはグラフを描画して見てみることはデータの特性を見るのにはいい手段だ。今回は株価データを単純に時系列プロットすることにする。

株価データを取得する

少し前だとyahoo financeやgoogle finance、あとは2017年で終わってしまったk-dbなどからデータを取ってきていたが今は現実的ではない。今だとQuandlというサイトからデータを取得するのがいいだろう。英語のサイトだがpythonから簡単なコードでデータを取得できるので使いやすい。


まずはライブラリのインポート

%matplotlib inline
import matplotlib.pyplot as plt
import quandl  #pipで別途インストールが必要
import pandas as pd

quandlがインストールされていない場合はターミナルで

>pip install quandl 

とやればインストールできる(おそらく)

株価データを取得する。今回は

の4つのデータを使うことにした。

データは次のコードで取得できる

df_quick = quandl.get('TSE/4318')  # Quick
df_hitachi = quandl.get('TSE/6501') # 日立製作所
df_takeda = quandl.get('TSE/4502')  #武田薬品工業
df_Nikkei = quandl.get('NIKKEI/INDEX')  # 日経平均株価

試しにQuickのデータを見ると

df_quick.head()
Date Open High Low Close Volume
2013-07-16 314.0 314.0 305.0 306.0 22100.0
2013-07-17 304.0 312.0 300.0 304.0 25700.0
2013-07-18 305.0 312.0 300.0 311.0 23300.0
2013-07-19 308.0 310.0 304.0 310.0 3600.0
2013-07-22 310.0 315.0 300.0 310.0 29100.0

株価をプロットする

のようになっている。今回は終値、つまり'Close'の値をプロットする。

# matplotlibの描画をかっこよく
plt.style.use('ggplot')

# 以下描画
fig = plt.figure(figsize=(18,8))
ax_quick = fig.add_subplot(221)
ax_hitachi = fig.add_subplot(222)
ax_takeda = fig.add_subplot(223)
ax_Nikkei = fig.add_subplot(224)

ax_quick.plot(df_quick.loc["2013-07-16": ,['Close']], label='Quick', color='b')
ax_hitachi.plot(df_hitachi.loc["2013-07-16": ,['Close']], label='HITACHI', color='r')
ax_takeda.plot(df_takeda.loc["2013-07-16": ,['Close']], label='Takeda', color='g')
ax_Nikkei.plot(df_Nikkei.loc["2013-07-16": ,['Close Price']], label='Nikkei', color='purple')

ax_quick.set_title('Quick')
ax_hitachi.set_title('HITACHI')
ax_takeda.set_title('TAKEDA')
ax_Nikkei.set_title('NIKKEI')

plt.show()

そしてこれが今回プロットした結果である。
f:id:doratai:20180202001816p:plain

どの銘柄も2016年の6月あたりから今に至るまで上がり続けている事が見て取れる。

ただ、縦軸のスケールが違うので、見た目の上がり方は同じでも、変動の幅は大きく違うことに注意したい。例えば、日立は2016年6月頃は400円ほどであり、2018年に入ることには900円程度になっているためその変動幅は約500円であるのに対し、武田薬品は2016年6月頃は約4000円であり、2018年に入る頃は6500円と、その幅は2500円ある。日立とは5倍の差がある。

銘柄ごとに株価の水準が異なるので、これをそのまま比較することはできない。そのためそのことを考慮した変動の割合を考える必要がある。




ただプロットしただけだけど今回はここまでにします。

学生時代にやっておくべきこと1位はなぜかいつも「海外旅行」

マ○ナビからメールが来た

社会人に聞いた!内定者が入社までにやっておいた方がいい遊びTOP10!
という記事がメールマガジンで紹介されてたので気になって見てみました。
そこには次のようにランキング形式でやっておいた方がいい遊びが紹介されていました。

1位 海外旅行をする 25人(9.6%)
2位 自分の好きな趣味に没頭する 17人(6.5%)
3位 国内旅行をする 13人(5%)
3位 合コン 13人(5%)
3位 飲み会 13人(5%)
6位 美術館・博物館めぐり 12人(4.6%)
7位 カラオケ 6人(2.3%)
7位 DVD鑑賞 6人(2.3%)
9位 映画鑑賞 5人(1.9%)
9位 アウトドア・スポーツ 5人(1.9%)
9位 キャンプ 5人(1.9%)

インドア派で休日は家から出ないのが基本な私にはなかなか難しい遊びが多い...

海外旅行       →そもそも日本から一回も出たこと無い。海外怖い。
趣味に没頭      →いつも思うんだけど趣味って何?
国内旅行       →お金無い。
合コン        →経験0。一度くらいは行ってみたいが呼ばれる可能性0。
飲み会        →どんだけ飲みたがりだよ...
美術館・博物館巡り  →これは面白そう。行ってみたい(博物館)。
カラオケ       →苦痛
DVD鑑賞       →寝る
アウトドア・スポーツ →インドアだっつてんだろ。
キャンプ       →あり得ない。

むしろ社会人になったほうがお金あるだろうし、旅行楽しそうだけどなあ...なんで学生時代に限定したがるんだ。
ふと、なぜこういうランキングでは海外旅行がいつも1位に来るんだろう?というどうでもいい疑問を抱いた。だってそこまでして行きたい理由分からないモン。必ず学生時代にやっておいたほうがいい事ランキングって必ず海外旅行1位だよね。もういいよ。そのランキングいらないよ。

理由1:「学生時代にやるべきこと=海外旅行」のイメージがつきすぎてる

いろいろなメディアで学生時代に海外旅行は必ず行くべき!みたいなことを言っている。これは先程のようなランキング形式のやつだったら絶対と言っていいほど出てくる。そういうものを普段から継続的に(この手の内容は毎年定期的に書かれてるよね?)目にすることで、多くの人に学生時代=海外旅行!みたいなイメージが既についてしまっている。

さて、アンケートの回答者は次のような質問を受ける。

「あなたが思う、学生時代にやっておいたほうがいいことはなんですか?」

ここで自分の周りの友達を思い出してほしい。そんなに皆が皆海外旅行に行ってる?少なくとも私の周りには片手で数えられるくらいしかいない。実際、学生の内に海外旅行に行く人達なんて一部だろう。その場合、アンケートの回答者は海外旅行に行ったことは無いが、「学生時代に海外旅行に行くべき!」という意識を刷り込まれた人達が必然的に多くなる。では先程の質問の回答はどうなるか?そう、

「海外旅行」

となるだろう。アンケートを回答する際、普通そんな真剣に考えない。そしたら当然こうなる用に思える。

理由2:社会人になって忙しくなる人達のほうが多い

当然だろ?と思うかもしれないがそうでないこともある。私は対して忙しく無いが、理系の大学に在籍しているため、特に実験などが多い研究室は週に休みがないのが当たり前みたいな教授の奴隷生活を送っている学生の話耳にすることもある。まあ、あくまで学生なので責任的な意味での大変さは無いだろうが時間的な意味での忙しさは、「社会人のほうが楽そう...」という言葉からなんとなく察する。まあ、ほとんど家から出ない私は低みの見物をするのみだが(紙とペンとパソコンと本あればなんとかなるのでそもそも研究室に行く意味ないし)。そんな私でもゼミのために勉強することは多いわ卒論書かなきゃだわと思ったより時間が無い(決して昔読んでた漫画やライトノベル等を読み返したり、アニメ見たりしているせいではない。そう、決して)。そんな中「時間のある学生時代の内に~」という言葉には違和感がある。学生時代はあくまで大学で「勉強」するための時間である以上海外「旅行」ではなく「留学」しろよ。どうせ海外いくなら。

大多数の人は社会人になってからのほうが時間的余裕がなくなるので、学生時代に海外旅行行っておけばよかったと後悔の念からそのような結果になるのだろう。でも思うんだよね。学生時代に海外行くようなアクティブな人なら、社会人になって忙しくても行くよね?

結局、学生時代も社会人になっても海外旅行行ってない人が答えている可能性...

理由3:1,2ときたら3がなきゃいけない気がする。

特に思いつかない。1、2ときたから一応様式美の3。

皆が学生時代に海外旅行に行くべきって言ってたから学生時代に海外旅行にいくべき。

つまりは、みんな(ソースは不明、ネットのランキングとか)が学生時代は海外旅行に行くべきって言ってたから学生時代は海外旅行にいくべきだと社会人が答え、それを聞いて社会に出た社会人がまた学生時代は海外旅行に行くべきって言ってたから学生時代は海外旅行にいくべきだと答え(以下無限ループ)により学生時代に海外旅行に行かなきゃならない感じになってる。

とりあえず学生時代にやるべきことランキングのまとめで

学生時代は海外旅行がおすすめ!!
海外旅行は学生時代のうちに!!

みたいなこと書くのはなんかなあ。旅行代理店か何かですか?それなら旅行プランへのリンク貼っとけよ。

注意

以上全て偏見でした。

カイ二乗分布の自由度を大きくしていくと同じ平均, 分散の正規分布と変わらなくなる話

カイ二乗分布

カイ二乗分布は大丈夫だと思いますが,  X_i \sim N(0, 1) のとき,

$$
\sum_{i=1}^n X_i^2 \sim \chi_n^2
$$

となります. 自由度 nカイ二乗分布に従うということです.

カイ二乗分布正規分布 N(n, 2n) に近づいていく様子

今回は自由度 n が大きくなるにつれて, カイ二乗分布が同平均, 同分散の正規分布に近づいていく様子をgifアニメ化してみました.

その結果がこちらです.

f:id:doratai:20171116135147g:plain

ここでは, 自由度を1から100まで動かしています. 自由度が大きくなるにつれて正規分布と一致していくのがわかります.

ソースコード

Jupyter 上で動かしました. アニメーションを作成するプログラムではなく, インタラクティブにnの値を変えることができ, そのときのカイ二乗分布(緑の破線)と正規分布(赤の曲線)の密度関数と, カイ二乗分布から得られた乱数をもとに作ったヒストグラムが描画されます.

%matplotlib inline
import matplotlib.pyplot as plt
from ipywidgets import interact
import numpy as np
from scipy.stats import norm, chi2

plt.style.use('ggplot')

bins = 100

@interact(n=(1,200))
def plot_chisquare(n=1):
    np.random.seed(0)
    data = np.random.chisquare(n, 10000)
    
    x = np.linspace(data.min(), data.max(), 1000) 
 
    param = norm.fit(data)
    pdf_fitted = norm.pdf(x, loc=param[0], scale=param[1])
    pdf = chi2.pdf(x,n)
    plt.plot(x, pdf_fitted, 'r-', label='normal pdf')
    plt.plot(x, pdf, 'g--', label='chi2 pdf')
    plt.hist(data, bins, color='pink', density=True)
    if n== 1:
        plt.ylim(0, 1)
    plt.title('Chi-square distribution with n degrees of freedom')
    plt.legend(loc='best')
    plt.show()

統計学を学ぶのにおすすめの問題集3冊

統計学を勉強しようと考えている方向けにレベル別におすすめの問題集を紹介する。統計検定やアクチュアリー等を考えている人にも十分なレベルを備えている問題集達なのでぜひ参考にしてください。

※統計検定準1級以上は範囲が広がるためここにある問題集だけでは少し不足かも...

[初級] 弱点克服 大学生の確率・統計 著者:藤田岳彦

統計学を勉強し始めたならまずはこの問題集に手を付けるのがいいと思います。高校生向けの問題集のようなデザインでとっつきやすく、確率統計の要点をしっかりおさえている、初学者にうってつけの問題集です。確率統計を勉強したことがある人なら必ず何処かで手にしたことがあるはず。初学者でなくとも、一度は通してやっておきたい一冊です。

[中級] 明快演習 数理統計 著者:小寺平治

弱点克服で確率統計に慣れたら次にやるべきはこの「明快演習 数理統計」です。難易度は弱点克服とそんなに変わりません。こちらのほうが、タイトル見れば分かる通り、数理統計が中心となっており、数理統計学の関する、いわば王道の問題を扱っています。受験的にいうと、数理統計学の頻出問題集のような感じです。統計学関連の問題を一通り網羅しているので、ここまでやればまず大学の数理統計学の授業で困ることは無いと思います。アクチュアリー、統計検定もここまでやればもう十分、といったレベル。

[上級] 確率統計演習2-統計- 著者:国沢清典

上記2冊をやってまだ物足りないという方がやるべきはこの一冊。結構レベルが上がります。数理統計の専門書すら扱ってないようなマニアックなものもあり、これ一冊でかなりの範囲を網羅可能。難しそうな見た目の割に全ての問題に解答がついており、問題集としてだけでなく参考書としても利用できるレベル。数理統計で難しめのレポート出されても、国沢統計を調べれば結構書けたりする。持っておくと結構役に立つ名著。ちなみにこれは「2」とありますが「確率統計演習1-確率-」もあります。

標本(不偏)分散の期待値, 分散[正規分布]

正規分布に従う確率変数の期待値, 分散等は統計に関連する本ならばまず間違いなく載っています. しかし, 標本(不偏)分散の期待値, 分散となってくるとなかなか取り扱っている本もサイトも少ない気がします. 定義から求めればいいといえばいいのですが, バカ正直に計算しようとすると結構大変です. 今回は標本分散の期待値, 分散について見ていこうと思います.

標本分散, 標本不偏分散の定義

標本分散S^2と標本不偏分散U^2を次のように書くことにします.

\begin{eqnarray}
S^2 &=& \frac{1}{n}\sum_{i=1}^{n} (X_i - \bar{X})^2 \\
U^2 &=& \frac{1}{n-1} \sum_{i=1}^{n} (X_i - \bar{X})^2
\end{eqnarray}

モーメント

期待値の計算に際し, モーメントを用いるので, 正規分布積率母関数と, 1次と2次モーメントを計算おこうと思います. 正規分布N(\mu, \sigma^2)に従う確率変数X積率母関数M_X(t) = \exp(\mu t + \sigma^2 t^2 /2)であるので,


E(X) = \frac{d}{dt}M_X(t)|_{t=0} = (\mu + \sigma^2 t) e^{\mu t + \sigma^2 t^2/ 2}|_{t = 0} = \mu
E(X^2) = \frac{d^2}{dt^2} M_X(t)|_{t=0} = \sigma^2 e^{\mu t+\sigma^2 t^2/2} + (\mu + \sigma^2 t)^2 e^{\mu t + \sigma^2 t^2/2}|_{t=0} = \mu^2 + \sigma^2

ちなみに, 標本平均\bar{X}\sim N(\mu, \sigma^2/n)より, そのモーメントは,
 \displaystyle
E(\bar{X}) = \mu \\
E(\bar{X}^2) = \mu^2 + \frac{\sigma^2}{n}

標本分散

まず初めに標本分散の方の期待値を見ていこうと思います.

期待値

期待値の定義から,

\begin{eqnarray}
E(S^2) &=& E \left[ \frac{1}{n}\sum_{i=1}^n (X_i - \bar{X})^2 \right] \\
&=& E\left[ \frac{1}{n} \sum_{i=1}^{n} X_i^2 - \bar{X}^2\right] \\
&=& \frac{1}{n} \sum_{i=1}^{n} E(X_i^2) - E(\bar{X}^2)
\end{eqnarray}

先程求めたモーメントを考慮すると,

 \displaystyle
E(S^2) = \frac{1}{n} \cdot n( \mu^2 + \sigma^2 ) - (\mu^2 + \frac{\sigma^2}{n})
       = \sigma^2 - \frac{\sigma^2}{n}
       = \frac{n-1}{n} \sigma^2
となる.

分散

続いて分散について見ていく. ここでは, \frac{\sum(X_i - \bar{X})^2}{\sigma^2} = \frac{nS^2}{\sigma^2} \sim \chi_{n-1}^2となる事実を用いる.
自由度n-1カイ二乗分布の分散が2(n-1)であることから,(環境の都合で \sigma^2 = s^2で表記します.)
\begin{eqnarray}
V(\frac{nS^2}{s^2}) &=& 2(n-1) \\
\frac{n^2}{s^4}V(S^2) &=& 2(n-1)\\
V(S^2) &=& \frac{2(n-1)s^4}{n^2}
\end{eqnarray}
よって(表記を戻すと)V(S^2) = \frac{2(n-1)\sigma^4}{n^2}になることが分かる.

※なんか,はてなでeqnarray環境使うとギリシャ文字がうまくいかないんだよなあ...

標本不偏分散

続いて標本不偏分散を見ていきます. といってもやることはほとんど同じ. 標本不偏分散のほうが期待値も分散もきれいになります.

期待値

先ほどと同様に定義からでも求まりますが, U^2 = \frac{n}{n-1}S^2の関係を用いれば,

 {\displaystyle
E(U^2) = E \left( \frac{n}{n-1} S^2 \right) = \frac{n}{n-1} E(S^2) = \frac{n}{n-1} \frac{n-1}{n} \sigma^2 = \sigma^2
}

標本不偏分散なので当然といえば当然.

分散

標本分散の時と同様に\frac{(n-1)U^2}{\sigma^2} \sim \chi_{(n-1)}^2であることを用いる.
{\displaystyle
V(\frac{(n-1)U^2}{\sigma^2}) = 2(n-1) \\
\frac{(n-1)^2}{\sigma^4} V(U^2) = 2(n-1) \\
V(U^2) = \frac{2\sigma^4}{n-1}
}
となる.

※eqnarray環境使わなかったため見た目が悪くなってます.


分散の分散を馬鹿正直にやると4次モーメントまで考えなくてはならなかったり, 式もかなり煩雑になるので, カイ二乗分布から持ってくる方が簡単で便利だと思います.
国沢統計を読んで見ると「4次モーメントを考えると...」という記述があってやろうとしたのですが途中で挫折しました...

Daftでグラフィカルモデルを作成してみる[Python]

森北出版の「Pythonで体験するベイズ推論」を読み進めていたら、2章で、Pythonのdaftというライブラリを用いて、グラフィカルモデルを作っていたのですが、そのソースコードは載っていなかったので自分で作ってみました。

作ったのは以下のグラフィカルモデル

f:id:doratai:20170605163828p:plain


参考にしたのは次のサイト。
daftでグラフィカルモデル
このサイトがかなり詳しく説明してくれています。

ソースコード

import daft
from matplotlib import rc
rc("font", family="Ricty", size=15)
rc("text", usetex="True")

pgm = daft.PGM(shape=[6,6])

# Nodes
pgm.add_node(daft.Node("alpha", r"$\alpha$", 4, 5)) # 名前 ラベル 座標
pgm.add_node(daft.Node("tau", r"$\tau$", 1, 4.5))
pgm.add_node(daft.Node("lambda_1", r"$\lambda_1$", 3, 4))
pgm.add_node(daft.Node("lambda_2", r"$\lambda_2$", 5, 4))
pgm.add_node(daft.Node("lambda",  r"$\lambda$", 2, 3))
pgm.add_node(daft.Node("obs", "obs", 2, 2, observed=True))

# Edges
pgm.add_edge("alpha", "lambda_1")
pgm.add_edge("alpha", "lambda_2")
pgm.add_edge("tau", "lambda")
pgm.add_edge("lambda_1", "lambda")
pgm.add_edge("lambda_2", "lambda")
pgm.add_edge("lambda", "obs")

pgm.render()
pgm.figure.savefig("pymc_p43.png")

多項分布

多項分布

基本性質

確率関数
 \displaystyle
f(x_1,\ldots,x_k) = \frac{n!}{x_1!\cdots x_k!} p_1^{x_1} \cdots p_k^{x_k}
期待値
 \displaystyle
E(X_i) = np_i
分散
 \displaystyle
V(X_i) = np_i(1-p_i)
共分散
 \displaystyle
Cov (X_i,X_j) = -np_i p_j

確率関数

1回の試行でk通りの可能な結果A_1,\ldots,A_kのいずれか1つのみが生じ、P(A_i) = p_i(i = 1,\ldots,k)とする。この試行を独立にn回繰り返したときに、A_iが生じる回数をX_iとするとき、X_1,\ldots,X_kの同時分布を多項分布といい、その確率関数は以下で与えられる。
$$
f(x_1,\ldots,x_k) = \frac{n!}{x_1!\cdots x_k!} p_1^{x_1} \cdots p_k^{x_k}
$$

期待値

\begin{eqnarray}
E[X_i] &=& \sum_{i=1}^k x_i \frac{n!}{x_1!\cdots x_k!} p_1^{x_1} \cdots p_k^{x_k} \\
&=& \sum_{i=1}^k \frac{n \cdot (n-1)!}{x_1!\cdots (x_i - 1)! \cdots x_k!}p_i \cdot p_1^{x_1} \cdots p_i^{x_i-1}\cdots p_k^{x_k} \\
&=& np_i \sum_{i=1}^k \frac{(n-1)!}{x_1!\cdots (x_i - 1)! \cdots x_k!} p_1^{x_1} \cdots p_i^{x_i-1}\cdots p_k^{x_k} \\
&=& np_i
\end{eqnarray}

分散

分散を求めるためにまず次の期待値を計算する。
\begin{eqnarray}
E[X_i(X_i - 1)] &=& \sum_{i=1}^k x_i(x_i - 1) \frac{n!}{x_1!\cdots x_k!} p_1^{x_1} \cdots p_k^{x_k} \\
&=& \sum_{i=1}^k \frac{n(n-1) \cdot (n-2)!}{x_1! \cdots (x_i-2)! \cdots x_k!} p_i^2 \cdot p_1^{x_1} \cdots p_i^{x_i-2} \cdots p_k^{x_k} \\
&=& n(n-1)p_i^2 \sum_{i=1}^k \frac{(n-2)!}{x_1! \cdots (x_i-2)! \cdots x_k!} p_i^2 \cdot p_1^{x_1} \cdots p_i^{x_i-2} \cdots p_k^{x_k} \\
&=& n(n-1)p_i^2
\end{eqnarray}
これと、V(X) = E(X(X-1)) + E(X) - E(X)^2であることを用いて分散を求める。
\begin{eqnarray}
V[X_i] &=& E[X_i(X_i-1)] + E[X_i] - E[X_i]^2 \\
&=& n(n-1)p_i^2 + np_i - n^2p_i^2 \\
&=& np_i(1-p_i)
\end{eqnarray}

共分散

期待値、分散のときと同様に計算することで、E(X_i X_j) = n(n-1)p_ip_jが得られるので、

 \displaystyle
\begin{eqnarray}
Cov(X_i, X_j) &=& E(X_i X_j) - E(X_i)E(X_j) \\
&=& n(n-1)p_i p_j - np_i \cdot np_j \\
&=& -np_i p_j
\end{eqnarray}

多項分布の例

サイコロをn回振ったときに、1の目がでる回数をX_1回。2の目がでる回数をX_2回。・・・6の目がでる回数をX_6回とする。また、それぞれの出る確率をp_i(i = 1,2,\ldots,6)とする。この時、確率ベクトル\boldsymbol{X} = (X_1,\ldots, X_6)'は、多項分布\mathrm{Multi}(n,\{p_i\})に従う。


参考文献

岩沢宏和(2012):『リスクを知るための確率・統計入門』,東京図書

「分かりやすい説明」という聞き手(読み手)の怠慢

「説明力」を求められる理系

最近、世の中から、特に理系に求められる能力として説明力がある。Wikipediaのような集合知、そこそこ専門的な知識でもググれば、その内容(理解できるかは別として)をすぐにでも確認することのできる時代に、専門用語をならべて偉そうに話す専門家の存在意義は薄れてきている(専門家が不要な訳ではない)。私自身、このことに対して、強く賛同している。分かりやすく説明できるということは、その分野に対して深い理解を持っていることの証明でもあるし、そのような説明ができる人間になりたいと考えてる。

 

氾濫する超入門書、内容の薄い○○でも分かる△△学

分かりやすさを求めるのは最早、世の中の大きな流れとなっていることが、書店に並ぶ本を見れば分かる。大きな書店に行けば、理工書がコーナーとしてあると思うが、そこに並ぶ本の中に少なからず、超入門系の本が置かれている。○○超入門!とか、文系でも分かる○○学!などといった本たちだ。目を引きやすい、いかにも優しそうなイメージの表紙に、「もう挫折しない!」みたいな帯がまかれていたりする。これらの本が一定数限りのあるならんでいるということは当然それなりに需要があるからならんでいるのだ。つまり、世の中がそのような本を求めているから、数少ない書店のスペースに、内容の薄さに対して、無駄に厚い超入門書が置かれるのである。

 

これらの本が、本当にわかりやすければいいのだ。数式を用いず、統計学ならば、各手法について、どうしてその手法なのか、どうやって使うのか、その結果は何を意味するのかを懇切丁寧に書いてくれるような書籍が本当の「分かりやすい」本なのではないかと思う。しかし、実際はどうだろう。特に私がよく見るのは統計関連のものだが、この手の本はだいたい、Excelの使い方を指南して終わる。なんか「回帰分析で予想!」みたいなことが書いてあって、Excelで一生懸命回帰分析の結果がでるまでの様子を、写真付き(これが、内容は薄く、本は厚くなる理由である)で説明するのだ。このことにどれほどの価値があるというのだろうか。正直、グーグルで検索したほうがよっぽどいい解説が出てくるだろう。そのレベルの無価値な本に、1600円、2000円なんて、どうして払えるだろうか?それ以上に、限りある書店スペースをそれらの本が占有していることが腹正しい限りだ。

「今」の入門書と「昔」の入門書の違い

分かりやすさが叫ばれるようになり、入門書のイメージも変わってきている。昔の入門書というのは「事前知識なくとも、読み進めることができる」というのが、「入門」の意味であった。例えば、多変量解析の入門書を読み進めるにあたって、行列の知識がかなり必要になってくる。昔の入門書というのは、このような必要な知識を付録として巻末に書いておいてくれる。また、当然、確率分布の定義から始まり、確率密度関数、期待値、分散共分散行列と、一つ一つ必要な知識を厳密に書かれているものであった。その分野の基本を詰めたようなものが入門書の意味するところであったのだ。それに対して、現在の入門書というのは、全く知識がない人に「なんとなくわかった気にさせる」ことが(超)入門書の意味するところになっている気がする。こういうものがあって、これはこう使います!といった、中身をそぎ落としすぎて、用語とExcelの使い方を説明する、「知ったか」を育成することを目的としているんじゃないかと疑いたくなるようなものだ。そもそも、学問というものはそれぞれの分野に一生をかける研究者がいる以上、本一冊読んだ程度で理解できるわけがない。

 

分かりやすさを求めた結果なくなるもの

分かりやすい説明をするためにはそのトレードオフとして「厳密性」が失われる。このことは多くの場合において重要ではないと考えられるが、それが失われることにより、多くの勘違いを生む結果となる。学問における定理、公式には、それが成り立つための「前提条件」というものが往々にして存在する。当然、説明する側はそのことは十分に理解しているだろうし、当然説明するときも、どこかしらで、その前提に触れていることだろう。しかし、聞き手側はどうだろうか。世の中の多くは、理系的な考え方に触れていない人が多く、そのような簡易な説明を求める層の多くはそのような人たちに占められると考えられる。そうなったときに、でかでかとプレゼンテーションで定理の内容が書かれていたのを見た聞き手が前提条件など覚えているだろうか。多くは覚えてなどいないだろう。関心があるのは、その結果にあるのだから。そうなってくると、当然前提条件が無視された定理の乱用が始まる。無価値な分析が世の中にはびこるようになる。統計等はその手法(Excelでなんかすれば結果出てくるし)よりも解釈のほうが難しい。そして、"Excelで分析した感を出した何か"を一生懸命発表するのだ。その結果が本当に正しいか、解釈があっているのかは二の次である。

 

何事も理解するためにはそれなりの努力が必要

分かりやすい説明を求めることは悪いことではない。ビジネスマンは研究者ではないので、少ない時間で概要をつかむ必要があるだろうし、研究者はそれにこたえられるほうがいいに決まっている。学生は、小難しい説明を永遠とされたらその授業に出る気など起きないし、勉強する意欲を失ってしまうだろう。

 

しかし、何を理解するにも、まっさらな状態で、なんとなく聞いていたのでは、いくらわかりやすい説明でも、理解できるわけがないのだ。その説明している人は数年から数十年をその分野にかけている人だろう。その内容の一部を短時間で、興味を持ってもらえるよう、もしくは、あなたが求めている部分だけに絞って説明してくれているのかもしれない。その努力に敬意を払い、自分なりに、理解する努力する必要がある。教えてもらうといった受け身思考ではよくない。学ぶんだという積極性が、たとえ、ガンガン質問にいくといった積極性までは持てなくとも、心持ちだけは積極でなくてはならない。そうでないと「分かりやすい」を追い求めた結果「内容の薄い」本や説明しか残らなくなるだろう(さすがにそんなことはないと思うが)し、そのような説明を聞いたり、その程度の本を読んだところで何も得るものなどない。そうならないためには、聞き手、読み手側にもそれなりの準備、努力が必要だと考える。

 

 

50年くらい前の本と現在本屋にならんでいる本を比べて、最近理論部分をしっかり扱う本が少ない(無いわけではない)気がして、ちょっとした愚痴でした。なんかすいません。

ベータ関数,ベータ分布

ベータ関数

定義

p,qが正の定数のとき、下記右辺の定積分を、p,qの関数と考え、ベータ関数と呼ぶ。

$$
B(p, q) = \int_0^1 x^{p-1} (1-x)^{q-1} dx \;\;\; (p,q > 0)
$$

ベータ関数とガンマ関数の関係

ベータ関数と、ガンマ関数の間には次の関係がある。
$$
B(p,q) = \frac{\Gamma (p) \Gamma (q)}{\Gamma (p + q)}
$$

ベータ分布

ベータ分布区間(0,1)上の確率分布であり、以下の確率密度関数によって定義される。

$$
f(x) = \left\{
\begin{array}{cc}
\frac{1}{B(p, q)} x^{p-1} (1-x)^{q-1} & (0 < x < 1) \\
0 & その他
\end{array}
\right.
$$
ベータ分布はBe(p, q)で表す。

以下X \sim Be(p,q)とする。

ベータ分布の平均,分散

\begin{eqnarray}
E[X] &=& \frac{p}{p + q} \\
V[X] &=& \frac{pq}{(p+q)^2 (p+q+1)}
\end{eqnarray}

導出

まずは平均について
\begin{eqnarray}
E[X] &=& \int_0^1 x \frac{1}{B(p,q)} x^{p-1} (1-x)^{q-1} dx \\
&=& \frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \int_0^1 x^p (1-x)^{q-1} dx \\
&=& \frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \frac{\Gamma(p+1)\Gamma(q)}{\Gamma(p+q+1)} \\
&=& \frac{\Gamma(p+q)}{\Gamma(p+q+1)} \frac{\Gamma(p+1)}{\Gamma(p)} \\
&=& \frac{p}{p+q}
\end{eqnarray}
分散を求めるに当たって、次のモーメントを求める。
\begin{eqnarray}
E[X(X-1)] &=& \int_0^1 x(x-1)\frac{1}{B(p,q)} x^{p-1} (1-x)^{q-1} dx \\
&=& -\frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \int_0^1 x^p (1-x)^q dx \\
&=& -\frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \frac{\Gamma(p+1)\Gamma(q+1)}{\Gamma(p+q+2)} \\
&=& -\frac{pq}{(p+q+1)(p+q)}
\end{eqnarray}
これより、V(X) = E(X(X-1)) + E(X) - E(X^2)から、分散は
\begin{eqnarray}
V[X] &=& E[X(X-1)] + E[X] - E[X]^2 \\
&=& -\frac{pq}{(p+q+1)(p+q)} + \frac{p}{p+q} - (\frac{p}{p+q})^2 \\
&=& \frac{pq}{(p+q)^2(p+q+1)}
\end{eqnarray}

モーメント

ベータ関数のk次モーメントを求める。
\begin{eqnarray}
E[X^k] &=& \int_0^1 x^k \frac{1}{B(p,q)}x^{p-1}(1-x)^{q-1} dx \\
&=& \frac{1}{B(p,q)} \int_0^1 x^{p+q-1}(1-x)^{q-1} dx \\
&=& \frac{B(p+k, q)}{B(p,q)} \\
&=& \frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \frac{\Gamma(p+k)\Gamma(q)}{\Gamma(p+q+k)} \\
&=& \frac{\Gamma(p+q)\Gamma(p+k)}{\Gamma(p)\Gamma(p+q+k)}
\end{eqnarray}

p,qが正の整数のときは

 \displaystyle
E [X^k] = \frac{(p+q-1)! (p+k-1)!}{(p-1)!(p+q+k-1)!}
と書くことができる。

ベータ分布の密度関数のグラフ

ベータ分布のグラフは、パラメータごとに以下のようになる。
f:id:doratai:20170528192802p:plain
pythonで描画。以下のサイトをおおいに参考にした。
【Python】scipyとmatplotlibでベータ関数を描画 - 歩いたら休め
ソースコードは以下。

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

x = np.linspace(0, 1, 1000)
plt.xlim(0,1)
plt.ylim(0,5)
plt.xlabel(r"$x$", fontsize=20, fontname='serif')
plt.ylabel(r"$f(x; p,q)$", fontsize=20,fontname='serif')
plt.title("PDF of Beta Distribution")
params = [[3,9],[6,6],[9,3],[1,1],[1,5],[5,1],[16,16]]
for param in params:
    rv = scipy.stats.beta(param[0], param[1])
    y = rv.pdf(x)
    plt.plot(x,y,'-',lw=2,label=param)
    plt.legend(bbox_to_anchor=(1.05,1), loc='best', borderaxespad=0) #凡例を枠外表示
plt.show()
ベータ関数の特徴

ベータ関数の特徴として、上記のグラフを見ればわかると思うが、パラメータp,qの値によっていろいろな形を取る、ということがある。例えば(p,q) = (1,1)のときは一様分布になっていることがグラフからもわかる。

ベータ分布とベイズ

ベータ分布がよく出てくるのは、ベイズ統計の分野である。ベイズ統計では、事前分布と、事後分布というものを考えるが、その時の事前分布としてベータ分布はよく使われる。グラフでも示した通り、ベータ分布のグラフはパラメータによって非常に柔軟に形を変えることができる。また、グラフとパラメータの対応関係をよく見て欲しい。(p,q) = (3,9)の時、グラフは左に偏り、(p,q) = (9,3)の時、グラフは右に偏っている。つまり、p,qの比がそのままグラフに表現されるのである。このことは、確率的主観を表現する際に都合が良い。予想が6-4であるとして、期待値を0.6とするならば、「私の確信」をp:q = 6:4のベータ分布で表現できるのである。ベイズ主義は確率を「ある事象をどれくらいできるか」の指標と解釈しているため、主観的な確率というものが非常に重要になってくる。その主観的確率を表現する際にベータ分布は非常に都合の良い分布なのだ。

参考文献

松原望,縄田和満,中井検裕(1991):『統計学入門(基礎統計学I)』,東京大学出版会
日本統計学会(2013):『日本統計学会公式認定 統計検定1級対応 統計学』,東京図書

一般逆行列の定義と存在

大学教養レベルで扱う線形代数では、逆行列は「正則行列(非特異行列)」である必要があり、\mathrm{rank}がフルランクであることが逆行列を持つ必要十分条件であった。しかし、行列が特異(逆行列を持たない、フルランクでない)である場合でも、逆行列を持つように、逆行列を拡張した、一般逆行列というものが存在する。統計学の中でも多変量解析などの分野では行列を多用するため、行列の話題というのは非常関心の高いものになる。そのような分野における、一般逆行列の利用は今や日常茶飯事らしいので、定義くらいは知っておきたい。

定義

m \times n行列\boldsymbol{A}の一般逆行列(generalized inverse)とは
$$
\boldsymbol{AGA} = \boldsymbol{A}
$$
を満たす、任意のn \times m行列\boldsymbol{G}のことである。
※一般逆行列の他に擬似逆行列、条件付き逆行列という用語で呼ばれることも多い
実際、\boldsymbol{A}が非特異である場合は\boldsymbol{G} = \boldsymbol{A}^{-1}であるので、\boldsymbol{AGA} = \boldsymbol{AA}^{-1}\boldsymbol{A} = \boldsymbol{A}となっている。

一般逆行列の存在

次に気になるのが、この「一般逆行列」が存在するのかどうかである。結論としては

あらゆる行列は少なくとも1つの一般逆行列を持つ。

これは次の定理で証明される。

定理

\boldsymbol{B}を最大列階数のm \times r行列。\boldsymbol{T}を最大行階数のr \times n行列とする。この時、\boldsymbol{B}は左逆行列\boldsymbol{L}を持ち、\boldsymbol{T}は右逆行列\boldsymbol{R}を持つ。そして、\boldsymbol{RL}\boldsymbol{BT}の一般逆行列である。

証明

\boldsymbol{B}が左逆行列\boldsymbol{L}を持ち、\boldsymbol{T}が右逆行列\boldsymbol{R}を持つことは既知としよう。この時、一般逆行列の定義から
$$
\boldsymbol{BT}(\boldsymbol{RL})\boldsymbol{BT} = \boldsymbol{B}(\boldsymbol{TR})(\boldsymbol{LB})\boldsymbol{T} = \boldsymbol{BT}
$$
である。すなわち、\boldsymbol{RL}\boldsymbol{BT}の一般逆行列である。
今、任意のm \times n行列\boldsymbol{A}を考える。\boldsymbol{A} = \boldsymbol{0}であるならば、明らかに任意のn \times m行列は\boldsymbol{A}の一般逆行列である。\boldsymbol{A} \neq \boldsymbol{0}であるならば、\boldsymbol{A} = \boldsymbol{BT}を満たす最大列階数の行列\boldsymbol{B}と最大列階数の行列\boldsymbol{T}が存在する。したがって、「あらゆる行列は少なくとも1つの一般逆行列を持つ」という結論を得る。


この証明で既知として用いている部分に関しては参考文献を参照ください。

参考文献

David A. Harville,(監訳)伊里正夫(2012) : 『統計のための行列代数(上)』,丸善出版

無限のパラドクス〜数学から見た無限論の系譜〜を読んで

理工系の新書として有名なレーベルの「BLUE BACKS」の本で、足立恒雄先生の著書、「無限のパラドクス」を読みました。

 

この本は、現代の無限論にたどり着くまでの歴史的経緯について非常に簡潔に分かりやすく書かれており、高校レベルの数学の知識があれば難なく読むことが出来ると思います。

 

無限論とは別に読んでいて感動したこと(ライプニッツの夢)

この本の中盤で、ライプニッツ(1646〜1716)の話が出てきます。この本によると、論理学の分野で「AならばB」のように、一般命題に文字をわりあてるのはライプニッツが始めたことだそうです。ライプニッツが夢見たのは「全ての論証を記号化し、算術問題に還元すること」。これを実現すると、裁判などで、データを投入するだけで後は計算機が計算して公平な判決を下してくれる。このような事を夢見ていたそうです。

少し前ならば、なんて夢物語なのだろうと私も思いました。しかし、人工知能、バックデータ解析といった用語がトレンドな現代、ライプニッツの夢見たことは現実味を帯びてきています。将来的にコンピュータに仕事を奪われると予想される職業に「弁護士」が入っているのですから。

 

人工知能が判決を下す。なんか、少し前のノイタミナのアニメで「PSYCHOPATH」という作品が人気を博しましたが、それを思い出します。警官の持つ拳銃型の「ドミネーター」という武器が、対象者の犯罪係数を算出し、その数値を元にその場で刑が執行されます。この犯罪係数が高いとその場で死刑執行も行われます。この世界はコンピュータ(最後にその中身も明かされる)に支配された人間の物語とも言えます。興味がありましたら見てみてください。映画化もされていた...かと思います(多分)。

 

科学が発達し、SF作品なども豊富な現代でこそ、想像のつくような世界を400年前の偉人が夢見ていたというのは想像もつきません。現代に生きる人々の中にも、私たちが理解出来ない様な未来をイメージして研究に取り組む天才もいるのかもしれないと思うと、私たちが夢物語だと思っていることも、500年も経てば現実化してるんでしょうか。最も最近の技術革新は10年前には映画の中の事だったことが当たり前のように行われるような進歩の仕方をしているので、それ以上となるとやはり想像できません(笑)。

 

なんだか、最終的に書評でもなんでもない所に帰着しましたが、興味深い本でしたので是非オススメしたいと思います。

 

また、著者の無限に関する本として別に「無限の果てになにがあるか」というものもありますのでそちらも読んでみてください。共に興味深い内容となっています。

日記的なの「統計の誤用を防ぐ書籍」

本日新宿の紀伊国屋に立ち寄りました。現在紀伊国屋では

 

「ダメな統計学を防ぐための書籍」

 

フェア?をやってるのか、理工書の階に行くと、中央のわかりやすい所にコーナーが設置されています。そこで見つけたフリーペーパーで「ダメな統計学を防ぐための書籍11冊」が紹介されていました(まあ、そのコーナーに置いてある本なんですが)どれも欲しい本ばかりなのですが、なにぶん学生故にお金が無いため買うことが出来ません...

 

あぁ、早く就職してお金が欲しい今日このごろです。

ちなみに、そのフリーペーパーの内容はブログに掲載されてる内容なので、次のリンクをクリックしてくれれば飛びますので(∩´。•ω•)⊃ドゾー

http://id.fnshr.info/2017/01/27/no-sdw/:ダメな統計学 11冊