はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その2
前置き
環境構築
まずPythonを実行する環境を作ります。 一応メモとして書きますが、他にもPythonの導入の記事はたくさんありますので、そちらも参考ください。 以下のものはOSXを想定しています。
自分は基本的にはpyenv + pyenv-virtualenvを利用しています。pyenvで切り替えつつ、学習する本の内容に合わせてpyenv-virtualenvで環境を作っています。 導入は各リポジトリのページに書いてあるとおりです。
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章 データの整理とベイズの定理
内容
一応箇条書きで書いていますが、メモです。 ちゃんと勉強したいときは、本文を参考にしてくださいね。
- 分布の紹介
- 正規分布
- 一様分布
- データの解釈
章末問題
- 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だとめんどくさかったりしますね。 特定の幅長での度数表や、最頻値について簡単に求めることが出来る方法を知っていらっしゃる方が居たら、是非コメントを頂けると助かります!