Easy to type

非正規労働者の職業訓練記録です。ボーナスと福利厚生を勝ち取る夢を持っています。

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その2

前置き

環境構築

まずPythonを実行する環境を作ります。 一応メモとして書きますが、他にもPythonの導入の記事はたくさんありますので、そちらも参考ください。 以下のものはOSXを想定しています。

自分は基本的にはpyenv + pyenv-virtualenvを利用しています。pyenvで切り替えつつ、学習する本の内容に合わせてpyenv-virtualenvで環境を作っています。 導入は各リポジトリのページに書いてあるとおりです。

github.com github.com

pythonは3系を使います。今回利用するライブラリはそんなに変なものではないので、2系である必要はありません。 解析はスクリプトでやるよりも、jupyterを使う方が僕は好きです。Stan言語をいじれない問題点はありますが、パラメータを変えるのが楽です。

pyenv install 3.6.1
pyenv virtualenv 3.6.1 3.6.1_postp
mkdir -p ~/Desktop/postp # 作業ディレクトリは適当に変えてください。
cd ~/Desktop postp
pyenv local 3.6.1_postp
pip install numpy scipy pandas matplotlib cython seaborn jupyter
pip install pystan

1章 データの整理とベイズの定理

内容

一応箇条書きで書いていますが、メモです。 ちゃんと勉強したいときは、本文を参考にしてくださいね。

  • 分布の紹介
  • データの解釈
    • なぜ正規分布にあてはめる?
      • 観測データの実際の分布は、より複雑な謎の分布に従うと予測される
      • 正規分布はパラメータが少なく、扱いが容易
         * データが分布にあてはめられることで、分布の数式を元にしたマニピュレーションが可能
        
    • どのようにあてはめる?
      • 頻度主義 : 標本平均、標本偏差を推定して、あてはめる
      • ベイズ統計 : データからパラメータの分布を推定し、その曖昧さを踏まえて分布を推定
        • ベイズの定理より、データが得られた時のパラメータ分布が得られる
                * ベイズの定理のために必要なもの
          
          • 尤度
          • 事前分布
            • 無情報事前分布の利用
          • 正則化定数(扱いが難しい!)
          • そこでMCMCですよ。

章末問題

  • pythonにて解析をしていきます。

データの準備

僕の好みですが、printではなくloggingモジュールを利用して途中経過を書いていきます。計算時間なども把握できるのが便利です。

import numpy as np
from logging import getLogger, Formatter, StreamHandler, DEBUG

# printではなくloggerを使う
def get_logger():
    logger = getLogger(__name__)
    logger.setLevel(DEBUG)
    log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s'
    formatter = Formatter(log_fmt)
    stream_handler = StreamHandler()
    stream_handler.setLevel(DEBUG)
    stream_handler.setFormatter(formatter)
    logger.addHandler(stream_handler)
    return logger

# データの準備
x = np.array([36,38,51,40,41,52,43,31,35,37,49,43,43,41,36,53,43,26,45,37,33,38,33,35,36,28,46,41,32,49,43,38,46,46,46,45,44,40,38,37,35,39,31,55,48,32,37,37,45,39,42,40,40,50,38,51,29,44,41,42,43,36,38,33,32,42,43,40,46,54,37,24,47,35,35,47,38,31,41,39,40,43,37,45,38,42,48,43,38,48,47,44,42,36,50,36,55,51,38,33])
logger = get_logger()

度数分布表

いきなりですが、pandasでは階級幅を決めるメソッドが無いので超めんどくさい… 手動でデータを入れるbinsを作成した後、pd.cutを使って仕分けます。

import pandas as pd
df = pd.DataFrame(x, columns=["weight"])

width = 5
minim_range = int((df.min() / width).apply(np.floor) * width)
maxim_range = int((df.max() / width).apply(np.ceil) * width)
bins = range(minim_range, maxim_range + 1, width)
c = pd.cut(df["weight"], bins)

d_df = df.groupby(c).count()
d_df.columns = ["count"]
d_df.columns.name = "range"
d_df.index.name = None 
d_df

ヒストグラムの作成

seabornを使って可視化します。matplotlibではないのは、やはり好みです pandas.plotと違ってseaborn / matplotlibのhistメソッドはbinsに配列を受け取ってくれます。

import seaborn as sns
import pandas as pd
%matplotlib inline
df = pd.DataFrame(x, columns=["weight"])

width = 5 
minim_range = int((df.min() / width).apply(np.floor) * width)
maxim_range = int((df.max() / width).apply(np.ceil) * width)
bins = range(minim_range, maxim_range + 1, width)
sns.distplot(df["weight"], kde=False, bins=bins)

統計量

最頻値だけはメソッドが無いので、countをとって最大の値を持つものを返しています。 pandasのvar, stdメソッドはデフォルトパラメータだと母分散、母標準偏差を返す事に注意してください。 numpyは標本分散、標本標準偏差を返す。謎の挙動。

import pandas as pd
df = pd.DataFrame(x, columns=["weight"])

# 標本平均
df_mean = df["weight"].mean()
logger.info("平均値は{0}".format(df_mean))

# 標本分散
df_var = df["weight"].var(ddof=0)
logger.info("標本分散は{0}".format(df_var))

# 標本標準偏差
df_std = df["weight"].std(ddof=1)
logger.info("標本標準偏差は{0}".format(df_std))

# データのソート
df["weight"].sort_values()

# 最大値
df_max = df["weight"].max()
logger.info("最大値は{0}".format(df_max))

# 最小値
df_min = df["weight"].min()
logger.info("最小値は{0}".format(df_min))

# 中央値
df_medium = df["weight"].median()
logger.info("中央値は{0}".format(df_min))

# 最頻値
df_most_freq = df["weight"].value_counts().iloc[0]
logger.info("最頻値は{0}".format(df_most_freq))

正規分部の扱い

scipy.statsのnormオブジェクトを利用します。 片側だけから特定の確率になる値を求めるメソッドが無いですが、intervalメソッドで両側の場合の値を求めて、2で割って考えれば良いです。 これは正規分布の対称性を利用したものですね。

import pandas as pd
from scipy.stats import norm 

df = pd.DataFrame(x, columns=["weight"])

# パラメータを設定
df_mean = df["weight"].mean()
df_std = df["weight"].std(ddof=1)
logger.info("正規分布のパラメータ: 平均は{0}、標準偏差は{1}".format(df_mean, df_std))

# 分布の平均が40.64なので、40付近の値が観測されやすい

# 45以上の値が観測される確率
# 1から45までの値が観測される確率を引く(1-cdf)
# scipyでは、sfを使えば良い
from scipy.stats import norm 
prob_upper_45 = norm.sf([45], loc=df_mean, scale=df_std)[0]
logger.info("45以上の値が観測される確率は{0}".format(prob_upper_45))

# 35以上45未満の値が観測される確率
# 45以下が観測される確率から、35以下が観測される確率を引く
prob_upper_35 = norm.cdf([35], loc=df_mean, scale=df_std)[0]
prob_upper_45 = norm.cdf([45], loc=df_mean, scale=df_std)[0]
prob_between_35_45 = prob_upper_45 - prob_upper_35
logger.info("35以上45未満の値が観測される確率は{0}".format(prob_between_35_45))

# 95%信頼区間
# 正規分布において、平均 ± 1.96 * 標準偏差の範囲が該当することを利用する
#min_edge = df_mean - df_std * 1.96
#max_edge = df_mean + df_std * 1.96
#prob_min = norm.cdf([min_edge], loc=df_mean, scale=df_std)[0]
#prob_max = norm.cdf([max_edge], loc=df_mean, scale=df_std)[0]
# prob_sig_range = prob_max - prob_min
# が、scipyだとintervalを利用すれば良い
min_edge, max_edge = norm.interval(0.95, loc=df_mean, scale=df_std)
logger.info("95%信頼区間は{0}から{1}".format(min_edge, max_edge))

# p(x>a) = 0.05であるような点a
# F(平均 + 1.64 * 標準偏差)が大体0.95になる
#prob_less_095 = norm.sf([df_mean + df_std * 1.64], loc=df_mean, scale=df_std)[0]
# 何も情報が無いとしたら、正規分布性を元に、90%信頼区間を出して、上端を利用すれば良い
_, less_095 = norm.interval(0.9, loc=df_mean, scale=df_std)
prob_less_095 = norm.cdf([less_095], loc=df_mean, scale=df_std)[0]
logger.info("a = {0} であれば、p(x>a)となる確率は{1}である".format(less_095, prob_less_095))

# 3つの4分位点
# 50%信頼区間と中央値を出せば良い
quad_first, quad_third = norm.interval(0.5, loc=df_mean, scale=df_std)
quad_second = df_medium
logger.info("第1四分位点は{0}, 第2四分位点は{1}, 第3四分位点は{2}".format(quad_first, quad_second, quad_third))

Rで出来ることがPythonだとめんどくさかったりしますね。 特定の幅長での度数表や、最頻値について簡単に求めることが出来る方法を知っていらっしゃる方が居たら、是非コメントを頂けると助かります!