統計・確率のお勉強

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

FACT FULNESS読みました

久しぶりの投稿です(いつも)。 2ヶ月くらい前?だったかに受けた応用情報処理技術者の合格通知がやっときました。一応IT企業に勤めているからと思って取得したのですが、普段コーディングとかしかしてないんであんまり内容が身近に感じないんですよね。業界的には標準の資格のような気がするのですが・・・。普段仕事で聞いたことあるような単語とかが出てきている(多分)のと、基本情報も春に取得したのもあって、そんなに時間かけず(6時間)に取得できたのでよかったです。筋トレとか初めてからいける気がする!って感覚が増えてきて、ちょい無謀気味な試験の受け方をしている気がする今日この頃です。まあ、結果がついてきてるんでいいですかね。FEもAPも大学生とかで取る人が多いんでそんな自慢する物でもないですが、他に話すところがないんだ、察しろ。

そんなことはどうでもいいとして、今日は最近読んだ本で面白かった本の紹介です。

FACT FULNESS

普段、書店行かれてる方なら多分見かけているんじゃないですかね? 少なくとも私がよく行く職場近くの書店では平積みされています。本のタイトルは最近流行りのマインド・フルネスからきています。マインドフルネスはApple Watchの機能としてあったり、メタ認知を高めたり、セルフコントロールを高めたり、近年注目されていますね。マインド・フルネスには人生の質を高めてくれる効果があります(ちゃんと研究もあるそうですよ)それと同じように私たちにはファクト・フルネスが必要だと、著者は説いています。

チンパンジーと私たち

著者は問いかけます。「チンパンジーと私たち、どちらが本当に世界を知っているかな?」と。ちなみに私はチンパンジーと同じくらい世界を知っていました。 「馬鹿にするな!チンパンジーより世界を知っているに決まっている!」と思った方は是非本書を手にとって冒頭の問題を解いてみてください。逆に、世界のことなんて全然分からないという謙虚な方は安心してください。自分のことを客観的にみれている証拠です。本を読み進めて、世界の事実知ってください。

ただ、この本は世界の事実を伝えることが目的の本ではありません。「事実に基づいて世界を見よう!」と著者が人生をかけて私たちに伝えるために作り上げた本です。

思い込みをやめて、データを見よう

私たちは自分が思っている以上に、経験と思い込みに支配されています。私の周りにはいませんが、「経験によると・・・」と言い出す人っていますよね?経験も大事です。それ以外に頼れない時もありますが、サンプル数1であることにたちかえれば、かなり怪しいですよね。自分の考えてることは本当に正しいのか、データに基づいて考える習慣をつけましょう。自分の論の根拠はどこにあるのかを考える癖をつけましょう。時間がかかるのでいつでもできるわけではありませんが、習慣づけることは大事です。著者もそれを望んでいます。だって世界的エリートたちですら世界のことを何も知らないのですから、知らなくても仕方ないと著者は語ります。その代わり、今から知ればいいのです。この本を読んで世界の事実を知り、データに基づいて考えるということを実感しましょう。

10の本能

著者は人が世界を間違ってみてしまう理由として以下の10個の本能をあげています。

  1. 分断本能
  2. ネガティブ本能
  3. 直線本能
  4. 恐怖本能
  5. 過大視本能
  6. パターン化本能
  7. 宿命本能
  8. 純化本能
  9. 犯人探し本能
  10. 焦り本能

それぞれがどういうことを言っているかは、実際に本を手にとってからのお楽しみです。

私がこの本を手にとった理由

私はこの本を発売前にAmazonで予約しました。普段から、月3万くらいは本につぎ込んでいるので、特別これがというわけではありませんが、この本に期待していたことは確かです。 私は大学で数理統計学を学んでいましたが、これをどう使うのかということがよく分かりませんでしたし、今もよくわかっていません。卒論では難しいことをやってみようと統計学の中では比較的難しい(主観)スパースモデリングに挑戦したりもしました。 しかし、やっぱり何のためにやってるのかよく分かりませんでした。この本にはデータをどう活用しているのかを学ぶつもりで購入しました。

もっとエビデンスに基づいた意思決定をする世の中になってほしいという思い

私は、あまり人と接するのが得意な方ではなく、交友関係も狭いので思い込みなのかもしれませんが、「それ、何に基づいて言ってるの?」って感じることがあります。テレビのコメンテーター、駅前で誰も聞いてないのに旗を掲げて、よくわからないことを話している政治家。最近はネットの力もあり、TVの影響力は薄れてきた感はありますが、日本全国に影響力を及ぼす場所で、思い込みや、間違った情報を垂れ流さないでほしいですし、国の行く末を決める政治家が感情論と人気取りで政策を決められては目も当てられません。それを防ぐためには、私たち一人一人が統計リテラシーを持つべきだと感じています。それ、根拠は?と突っ込むようになるだけでもだいぶ変わってくると思います。

この本は、上記の思いに通じることが多分に含まれると考えます。読んでみて、事実に基づいて考えるということ知っていただきたいです。

難しいことより、伝えること

この本の著者とその家族は、人に世界の真実を知ってもらおうと色々な工夫をしています。グラフィカルなバブルチャートについて、よくこの本の中で言及されます。 近年のAIブームもあり、最近は頭のいい人たちがkaggleとかで機械学習を競っているようです。私はやっていないのでよく分かりませんが。 あれはあれで応用が色々あるのだと思いますが、データ分析とはそんな難しいことをする必要はないんじゃないかと思います。だって、難しすぎて、伝わらないし。それに今はやってるDeep Learningは中身基本ブラックボックスですから。

そんなことよりも、簡単に作ったヒストグラムや折れ線グラフを見てなぜを考察する方が人にも伝えやすいし、有意義なんじゃないかなあと思います。

著者も多くの人に伝えるためにかなり骨を折ったようです。執念がすごい。

ビジネス書には珍しい筆者の思いが伝わってくる内容

この本は、世界中の人々に著者が事実に基づいて世界を見てほしいと願って作り上げた人生最後の作品です。「おわりに」の部分の内容を読み終わった時には、小説を読み終わった後のような虚無感があります。本の厚さに対して値段が少し安く設定されているような気がするのですが、これは筆者の多くの人に手にとってほしいという願いがあるのではないかと思ったりしたりもしました。本当に著者と、その家族の人生が詰まった本です。読んで後悔しません。ここまでいい本はなかなか巡り会えないと思います。

もし気になった方は、一度手にとっていただいたけたらと思います。

はてなブログでjavascriptが書けることを知った

足し算するだけのコードを試しに書いてみた。

最近仕事でjavascriptを書く機会があったこともあり、ふと思った。普段書いてるブログってjavascript埋め込めたら結構色々なことできね?と。 調べたらはてなブログjavascriptを埋め込むことができるらしい。 今回は実際に埋め込めるか試してみることにした。 今回はjqueryを使うので次のタグを最初に埋め込んでいる。

<script src="https://ajax.googleapis.com/ajax/libs/jquery/3.3.1/jquery.min.js"></script>

次に、入力とかするためのinputタグ(html)を使って入力できるようしたものである。

+ =

上記のを表示するためのhtmlは以下

<input id='plus1' type=number> + <input id='plus2' type=number > = <span id='result_plus'></span><br>
<input id='sum_button' type=button value='計算'>

数値を入力して計算ボタンを押せば、足し算された結果が表示される。(the サンプルプログラム)

※試しに書いただけのものなので入力チェックとかはしてません。

$(function() {
    let add = function(a, b) {
        return Number(a) + Number(b);
    }; 
    $('#sum_button').click(function() {
         a = $('#plus1').val();
         b = $('#plus2').val();
         res = add(a, b)
        $('#result_plus').text(res);
    });
});

jqueryで入力された値を取得して、ボタンが押されたら足してそれを対象箇所に出力している。

やってみた感想

はてなブログは当たり前だがブログ用なのでプログラムを書くようには記事編集のところはできていない。 今回は簡単なものだったので直接打ち込んだが、プログラミングエディタではないので非常に面倒臭かった。 はてなブログに埋め込む用のスクリプトは別のところで書いて、出来上がったものを貼り付けるようにすれば楽だし色々できそうだ。 Tensorflow.jsなんてものも出てきてるし、無料ブログで機械学習みたいなこともできそうだ。

気が向いた時の知識のアウトプットだったはてなブログの楽しみが一つ増えた。

ベイズ統計:パラメータの点推定

久しぶりなので少し短めに・・・
大学卒業して少し統計離れてましたけど、土日やることもないので、前々から気になっていたベイズ統計の勉強でも少しづつ始めていこうと思います。

パラメータの点推定

点推定
未知であるパラメータの値をデータから推定すること

通常の統計学なら不偏推定量とか、最小分散不偏推定量とか、十分統計量とかありました。それのベイズ版です。

ベイズ分析におけるパラメータの点推定は、パラメータに真の値と点推定の乖離度をある尺度(損失関数)ではかり、この損失関数が出来る限り小さくなるようにすることを考えます。

ここからは、パラメータ  \pi の点推定を \delta 損失関数を  L(\pi, \delta) と書くことにする。ベイズ分析によく使われる損失関数としては次のような関数があります。

\begin{align}
2乗誤差損失({\rm quadratic \, loss}) &: L(\pi, \delta) = (\pi - \delta)^2 \\
絶対誤差損失({\rm absolute \, loss}) &: L(\pi, \delta) = |\pi - \delta|^2 \\
0 - 1 損失 ({\rm 0-1 \, loss}) &: L(\pi, \delta) = 1 - \boldsymbol{1}_{\pi} (\delta)
\end{align}

ここで、  \boldsymbol{1}_{\pi} (\delta) \delta = \pi の時 1 それ以外のとき 0となる指示関数。

パラメータの真の値はわかっていないので、ベイズ分析の点推定では、損失関数  L(\pi, \delta) の期待値を  \pi の事後分布で評価したもの

\begin{equation}
R(\delta | D) = E_{p (\pi | D)} [L(\pi, \delta)] = \int_0^1 L(\pi, \delta) p(\pi | D) d\pi
\end{equation}

を考え、これを出来るだけ小さくするように点推定を選択します。 R(\delta | D)期待損失 (expected loss) と呼びます。これを最小化問題として定式化したものが、未知のパラメータ  \pi の点推定  \delta^*

\begin{equation}
\delta^* = \arg \min_{0 < \delta < 1} R(\delta | D) = \arg \min_{0 < \delta < 1} \int_0^1 L(\pi, \delta) p (\pi | D) d\pi
\end{equation}

となります。"  \arg " は 「  \delta^* \min_{0 < \delta < 1} R(\delta | D) という最小化問題の解である」という意味です。

参考文献

中妻照雄(2007):『ファイナンスライブラリー10 入門ベイズ統計学』, 朝倉書店

正規性の検定(Shapiro-Wilk検定)

正規性の検定

データを分析するにあたり、

  1. データが正規分布に従う
  2. データが独立な標本である

といった仮定を置くことは多い。そのような場合に分析をする際、これら二つの仮定が満たされているか確認する必要が出てくる。そのための手法として統計的仮説検定がある。今回はその中の Shapiro-Wilk検定Pythonで株価収益率に対してやってみようと思う。

二種類の誤り

仮説検定には2つの誤りがある。

第1種の誤り
帰無仮説が正しいにも関わらず、帰無仮説を棄却してしまう誤り。

第2種の誤り
帰無仮説が正しくないにも関わらず、帰無仮説を採択してしまう誤り。

通常、第一種の誤り確率をある水準(有意水準)以下に抑えた状態で、第二種の誤り確率をできるだけ小さくするようにする。そのため、帰無仮説が採択された場合は、積極的に帰無仮説を正しいとみているわけではなく、「帰無仮説が間違っているとは言えない」というような消極的な解釈となる。

Shapiro-Wilk検定

私が所持している本の中にはこの検定に関する記述のあるものが見つからなかったため、今回は
シャピロ–ウィルク検定 - Wikipedia
の内容でお茶を濁すことにする。

統計学における、シャピロ–ウィルク検定(シャピロ–ウィルクけんてい)とは、 標本  x_1, \ldots, x_nが正規母集団からサンプリングされたものであるという帰無仮説を検定する検定である。この検定方法は、サミュエル・シャピロとマーティン・ウィルクによって、1965年に発表された。検定統計量は、

$$
W = \frac{\left(\sum_{i=1}^n a_i x_{(i)} \right)^2}{\sum_{I=1}^n (x_i - \bar{x})^2}
$$

ただし

  • x_{(i)} は、i番目の順序統計量。
  • \bar{x} は標本平均。
  • 定数  a_i は次の式によって与えられる。

$$
(a_1, \ldots, a_n) = \frac{m' V^{-1} }{\left( m' V^{-1} V^{-1} m \right)^{1/2}}
$$

ただし、

$$
m = (m_1, \ldots, m_n)'
$$

m_1, \ldots, m_nは、標準正規分布からサンプリングされた独立同分布の確率変数の順序統計量の期待値であり、「V」は、この順序統計量の分散共分散行列である。帰無仮説は、「W」が小さすぎる場合に棄却される。

※一応Wikiの引用だが、少し書き換えてある。

兎にも角にもこれを用いて正規性を検定する。データは株式の収益率を使う。

Pythonでは次のようにすれば、Shapiro-Wilkの検定が可能である.

from spicy import stats
stats.shapiro(df['QUICK'])  # df['QUICK'] はデータフレームdfに入っているQUICKの収益率データを表す

結果は以下のようになった。

統計量W p値
QUICK 0.9260867834091187 1.0152270032202054e-22
日立 0.9507756233215332 1.2722628221227775e-18
武田薬品 0.9450637698173523 1.0990565268293669e-19

どれもp値がほとんど0に近いので、棄却される。つまり、どれも正規分布に従っているとは言えない事がわかる。

各業種と日経平均株価についてヒストグラムとガウシアンフィッティングで得られた曲線(もし、そのデータが正規分布に従うならばフィットするであろう正規分布な曲線)を一緒に描画した図をみてみよう。

f:id:doratai:20180207230849p:plain

これをみて見ると、正規分布よりも尖った分布になっているのがみて取れる。

参考文献の本では、4銘柄中3銘柄帰無仮説が採択されてたのに...

正規分布に従っているという仮定には無理がある事がわかった。

参考文献

[1] 横内大介, 青木義充(2014):『現場ですぐ使える時系列データ分析〜データサイエンティストのための基礎知識〜』技術評論社

収益率の相関を見る

相関係数

二つの銘柄X, Yの収益率データを \{x_1,\ldots, x_T \}, \{y_1,\ldots, y_T\} とする。この時相関係数

$$
corr(X,Y) = \frac{ \sum_{t=1}^T (x_t - \bar{x})(y_t - \bar{y}) }{\sqrt{\sum_{t=1}^T (x_t - \bar{x})^2} \sqrt{\sum_{t=1}^T (y_t - \bar{y})^2} }
$$

で定義される。今回は散布図行列をプロットして見ることにする。
今回は

の3銘柄で試してみた。

%matplotlib inline
import matplotlib.pyplot as plt
import quandl
import pandas as pd
import seaborn as sns
import numpy as np
from scipy.stats import norm

# データの取得

df_quick = quandl.get('TSE/4318')  # Quick
df_hitachi = quandl.get('TSE/6501') # 日立製作所
df_takeda = quandl.get('TSE/4502')  #武田薬品工業

# 対数差収益率に変換
Quick_RoR = np.diff(np.log(df_quick.loc["2013-07-16":, ['Close']]), axis=0)*100
Hitachi_RoR = np.diff(np.log(df_hitachi.loc["2013-07-16":, ['Close']]), axis=0)*100
Takeda_RoR = np.diff(np.log(df_takeda.loc["2013-07-16":, ['Close']]), axis=0)*100
Nikkei_RoR = np.diff(np.log(df_Nikkei.loc["2013-07-16":, ['Close Price']]), axis=0)*100

# データフレームに収益率をまとめる
df = pd.DataFrame(np.concatenate((Quick_RoR, Hitachi_RoR, Takeda_RoR), axis=1), columns=['QUICK', 'HITACHI', 'TAKEDA'])

# seaboard で散布図行列を出力
sns.pairplot(df, kind='reg',size=8)

得られた結果が以下

f:id:doratai:20180207003104p:plain

これを見る限り、武田薬品と日立の相関が比較的強いように見える。実際に相関を計算すると

df.corr()
QUICK HITACHI TAKEDA
QUICK 1.000000 0.324614 0.307317
HITACHI 0.324614 1.000000 0.511227
TAKEDA 0.307317 0.511227 1.000000

となる。日立と武田薬品は約0.5と相関があり、日立とQUICKと比べて大きい事がわかる。



少し短いけどここまでにします。

<時系列データ分析>ボラティリティを見る〜Pythonと株価データを使ってお勉強〜

ボラティリティ

ちゃんと本を読むまでボラティリティ のことを単純に標準偏差だと思っていたが、どうやら定義がいくつかあるようだ。詳しくはないのでWikipediaから引用する。

金融工学においてボラティリティ(volatility)とは、広義には資産価格の変動の激しさを表すパラメータ。

ボラティリティ - Wikipedia

細かいところはWikiをみてください。

とりあえずの理解としては株価のバラツキ具合だと思うことにする。株価のバラツキが大きいということは投資した場合に株価の値上がりで儲ける可能性も高いが、逆に失う可能性も高い。つまりはリスクが大きいということになる。ボラティリティは銘柄の特徴を捉えるのに重要な指標の一つのようだ。

ヒストリカル・ボラティリティ

ここではヒストリカル・ボラティリティの定義を確認する。ヒストリカル・ボラティリティは価格の対数差収益率の標準偏差で定義される。

t 時点の株価を P_t とした時の t 時点での対数差収益率を r_t とする。

$$
r_t = \log P_t - \log P_{t-1} = \log \frac{P_t}{P_{t-1}}
$$

そのため、ヒストリカル・ボラティリティn 個の収益率データ \{r_1, r_2, \ldots, r_n \} が得られているとき

$$
s = \sqrt{\frac{1}{n-1}\sum_{i=1}^n (r_i - \bar{r})^2 }
$$

で推定される。

対数差収益率をプロットする

各銘柄の変動のしやすさを視覚的に見るために対数差収益率をプロットしてみよう。
今回も前回同様銘柄は

の4つで行く。

import matplotlib.pyplot as plt
import quandl
import numpy as np

# データをquandlから取得
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')  # 日経平均株価

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

fig1 = plt.figure(figsize=(18, 10))
ax1_quick = fig1.add_subplot(221)
ax1_hitachi = fig1.add_subplot(222)
ax1_takeda = fig1.add_subplot(223)
ax1_Nikkei = fig1.add_subplot(224)

ax1_quick.plot(np.diff(np.log(df_quick.loc["2013-07-16":, ['Close']]), axis=0)*100, label='Quick', color='b')
ax1_hitachi.plot(np.diff(np.log(df_hitachi.loc["2013-07-16":, ['Close']]), axis=0)*100,label='HITACHI', color='r')
ax1_takeda.plot(np.diff(np.log(df_takeda.loc["2013-07-16":, ['Close']]), axis=0)*100 , label='Takeda', color='g')
ax1_Nikkei.plot(np.diff(np.log(df_Nikkei.loc["2013-07-16":, ['Close Price']])*100 , axis=0), label='Nikkei', color='purple')

ax1_quick.set_title('Quick')
ax1_quick.set_ylim=(-20,20)

ax1_hitachi.set_title('HITACHI')
ax1_hitachi.set_ylim(-20,20)

ax1_takeda.set_title('TAKEDA')
ax1_takeda.set_ylim(-20, 20)

ax1_Nikkei.set_title('NIKKEI')
ax1_Nikkei.set_ylim(-20, 20)

plt.show()

y軸の幅はQuickの変動幅が一番大きかったのを何回か実行する中で確認したので
それに合わせた。y軸の値は パーセント(%)

f:id:doratai:20180202170051p:plain

これを見るとボラティリティは Quick > 日立 > 武田薬品日経平均株価 となっている事が視覚的にわかる。

Quickは \pm 20 %と大きく変動しているが、日経平均武田薬品\pm 5 % の幅に収まっている。

つまり、この4つの中だとQuickが最もボラタイル(ボラティリティが大きい)である事がわかる。

データの要約

視覚的に見ることはできた。続いては、平均や分散のような統計量を求めて見ることにする。

収益率の平均

各銘柄の収益率の平均を求めてみよう。
次のようにすれば計算できる。

quick_mean = np.mean(np.diff(np.log(df_quick.loc["2013-07-16":, ['Close']]), axis=0)*100)
hitachi_mean = np.mean(np.diff(np.log(df_hitachi.loc["2013-07-16":, ['Close']]), axis=0)*100)
takeda_mean = np.mean(np.diff(np.log(df_takeda.loc["2013-07-16":, ['Close']]), axis=0)*100)
Nikkei_mean = np.mean(np.diff(np.log(df_Nikkei.loc["2013-07-16":, ['Close Price']]), axis=0)*100)

print("Quick : " + str(quick_mean))
print("HITACHI : " + str(hitachi_mean))
print("TAKEDA : " + str(takeda_mean))
print("Nikkei : " + str(Nikkei_mean))

結果は

銘柄名 QUICK 日立 武田薬品 日経平均
平均(%) 0.163 0.027 0.029 0.043

いづれの銘柄の平均値も正の値をとるが、どれも0に近く、小さな値になっている。

※計算結果(平均)
Quick : 0.162560953612
HITACHI : 0.0269943597075
TAKEDA : 0.0285518660969
Nikkei : 0.0426793473652

ボラティリティを計算する

続いては、各銘柄のボラティリティ標準偏差)を求める。
次のようにすれば計算できる。

quick_std = np.std(np.diff(np.log(df_quick.loc["2013-07-16":, ['Close']]), axis=0)*100)
hitachi_std = np.std(np.diff(np.log(df_hitachi.loc["2013-07-16":, ['Close']]), axis=0)*100)
takeda_std = np.std(np.diff(np.log(df_takeda.loc["2013-07-16":, ['Close']]), axis=0)*100)
Nikkei_std = np.std(np.diff(np.log(df_Nikkei.loc["2013-07-16":, ['Close Price']]), axis=0)*100)

print("Quick : " + str(quick_std))
print("HITACHI : " + str(hitachi_std))
print("TAKEDA : " + str(takeda_std))
print("Nikkei : " + str(Nikkei_std))

結果は

銘柄名 QUICK 日立 武田薬品 日経平均
ボラティリティ(%) 2.713 1.896 1.342 1.316

先ほどの図の印象と数字が一致している事がわかります。Quickのボラティリティは2.713% であり、武田薬品が 1.342% であることから大きさとしてはQuickは武田薬品の約2倍であり、より価格変動が大きい銘柄である事がわかる。


※計算結果(ボラティリティ)
Quick : 2.7130695301
HITACHI : 1.8958143491
TAKEDA : 1.34199583079
Nikkei : 1.31599308226

豆○ば〜

QUICKって、これを書く際に参考にしている「現場ですぐ使える時系列データ分析ーデータサイエンティストのための基礎知識ー」の著者が所属している会社なんだよ〜

<時系列データ解析> 株価をプロット〜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")