僕が統計を勉強するのに使ってる本まとめ(14/12/06更新)
統計を勉強し始めて3年ぐらいになりましたが、その3年間で「あー、これ何度も見直してるな」的な本をまとめておこうと思い立ったのが午後4時ぐらいなのですが、せっかくの機会なので、書いてみようと思います。ただ、とにかく本を読むのが苦手な僕なので、レパートリーに限りがありまして「偏ってる」のは間違いありません。
それから、数理科学をやってる人の本棚です（それ数学的にどうなん？とか気になっちゃうような人が書いてます）。数式見て、無理そうですという人は、あまり参考にされない方が良いかと思います。そういう方向けには、いろんな方がオススメの本のブログなどを書いているのでそちらを参考にされるのがいいかと。。。※あと、以下のAmazonリンクから購入していただいても、僕には1円も入りません（笑）。
数理統計学全般
Statistical Inference
- 作者: G.C. Casella,Roger L. Berger
- 出版社/メーカー: Brooks/Cole
- 発売日: 2008/06/07
- メディア: ペーパーバック
- 購入: 1人 クリック: 1回
- この商品を含むブログを見る
とっても重宝する統計学の教科書的立ち位置なのがこの本です。統計学で扱う「集合」の話から始めて、「分布」とは何か、「統計的検定」などなど、たくさん載っていて、数理統計の基本でつまづいたときに必ず見直す1冊です。
入門・演習 数理統計
- 作者: 野田一雄,宮岡悦良
- 出版社/メーカー: 共立出版
- 発売日: 1990/05
- メディア: 単行本
- クリック: 4回
- この商品を含むブログ (2件) を見る
こちらは和書です。上のものよりも、ちょっと分量が少ないですが、こちらも数理統計を勉強する上で重宝する一冊。
統計学のための数学入門30講
- 作者: 永田靖
- 出版社/メーカー: 朝倉書店
- 発売日: 2005/04/01
- メディア: 単行本
- 購入: 23人 クリック: 398回
- この商品を含むブログ (24件) を見る
統計学を勉強する前に、抑えておきたい数学の知識を整理している1冊です。といいつつ、僕も基礎が抜けててよくこの本に立ち戻ります。特に、行列計算などの基本的な部分や、逆行列の存在条件など、案外忘れてしまった知識を取り戻すのに役立ちます。
bayesian inference in statistical Analysis
Bayesian Inference in Statistical Analysis (Wiley Classics Library)
- 作者: George E. P. Box,George C. Tiao
- 出版社/メーカー: WILEY
- 発売日: 1992/03
- メディア: ペーパーバック
- この商品を含むブログを見る
ベイジアンモデリングなどをしていると、あまり数式が出てこなかったりして、計算が複雑なので結果だけ示されることがありますが、どうしても気持ち悪いとか、自分で導出したいぞ！という人向けの1冊です。事前分布・事後分布などの概念をそこそこ知ってないと手を出せないと思いますが、読めるようになれば楽しい1冊なんじゃないかなと。僕は好きです。
モデリング系
多変量解析概論
- 作者: 塩谷実
- 出版社/メーカー: 朝倉書店
- 発売日: 1990/03
- メディア: ハードカバー
- この商品を含むブログを見る
1990年に出版された本ですが、多変量解析の基本的な部分に立ち戻るのに重宝している1冊です。線形回帰モデルとは何かを数式と図を用いて丁寧に説明してくれますし、決定係数とは何かなど基本的なものをきちんと定義してくれています。
一般化線形モデル入門
- 作者: Annette J.Dobson,田中豊,森川敏彦,山中竹春,冨田誠
- 出版社/メーカー: 共立出版
- 発売日: 2008/09/08
- メディア: 単行本
- 購入: 15人 クリック: 152回
- この商品を含むブログ (13件) を見る
一般化線形モデルについて書かれた比較的読みやすい1冊かなと思っています。一般化線形モデルについては、有名なものとして久保先生の書かれている「データ解析のための統計モデリング入門」もあります（勉強会も開かれているほど、有名な1冊です）。ただ、僕は緑本は、もうちょっと数学的に書いてほしい。。。ということもあり、こちらの方をオススメします。
Generalized Linear Model, Second Edition
- 作者: P. McCullagh,John A. Nelder
- 出版社/メーカー: Chapman and Hall/CRC
- 発売日: 1989/08/01
- メディア: ハードカバー
- クリック: 47回
- この商品を含むブログ (4件) を見る
一般化線形モデルを提案したNelderが共著者になっている本です。僕もまだこの本全体を理解できていませんが、一般化線形モデルをかなり深くまで掘り下げて書かれている1冊です。ただ、結構わかりにくい気はしています。
Generalized Additive Models: An Introduction with R
- 作者: Simon Wood
- 出版社/メーカー: Chapman and Hall/CRC
- 発売日: 2006/02/27
- メディア: ハードカバー
- クリック: 17回
- この商品を含むブログ (8件) を見る
一般化加法モデル（一般化線形モデルに対して、スプラインを使って拡張する）についての本です。ここまでくると、かなり柔軟なモデリングが可能になります。一般化線形モデルはある程度わかっていないと読めないかもしれません。。。過去に、野球データ解析の際に用いています（下のスライドの後半を参照していただければ）。
20140727_第1回スポーツデータアナリティクス基礎講座
機械学習
The Elements of Statistical Learning
- 作者: Trevor Hastie,Robert Tibshirani,Jerome Friedman
- 出版社/メーカー: Springer
- 発売日: 2009/03
- メディア: ハードカバー
- 購入: 1人 クリック: 222回
- この商品を含むブログ (14件) を見る
入門パターン認識と機械学習
- 作者: 後藤正幸,小林学
- 出版社/メーカー: コロナ社
- 発売日: 2014/03/30
- メディア: 単行本
- この商品を含むブログ (5件) を見る
機械学習を数学やってる人が入門するのに結構いい本だと思いました。PRMLは個人的に読みにくいなーと思っていましたが、こっちは読みやすかったです。
Python, Rで分析・解析
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (6件) を見る
Pythonでデータ処理のやり方を覚えるのにまず最初にやった1冊です。この本で、大体の基本的なことを学びました。今でも見直して、「あー、そうだったそうだった」と想い出すのに使ってます。
集合知プログラミング
- 作者: Toby Segaran,當山仁健,鴨澤眞夫
- 出版社/メーカー: オライリージャパン
- 発売日: 2008/07/25
- メディア: 大型本
- 購入: 91人 クリック: 2,220回
- この商品を含むブログ (277件) を見る
Pythonで機械学習するなら、持っておいた方がいいかなと思う1冊です。とにかく、話題が豊富で、機械学習に興味のある人なら、どこかの章がひっかかるんじゃないかなと。
入門機械学習
- 作者: Drew Conway,John Myles White,萩原正人,奥野陽,水野貴明,木下哲也
- 出版社/メーカー: オライリージャパン
- 発売日: 2012/12/22
- メディア: 大型本
- 購入: 2人 クリック: 41回
- この商品を含むブログ (11件) を見る
とりあえず、機械学習をRでやってみたいというのであれば、この本が参考になりました。ただ、Rで学ぶデータサイエンスシリーズというのもありまして、こちらの方が好きな方もいらっしゃるようです。
Rで学ぶデータサイエンス（シリーズものです20巻ぐらいあるんじゃないかな）
- 作者: 藤井良宜,金明哲
- 出版社/メーカー: 共立出版
- 発売日: 2010/04/22
- メディア: 単行本
- クリック: 13回
- この商品を含むブログ (3件) を見る
こちらはシリーズ物で、Rで様々なデータ解析をやってみようというコンセプトのもとつくられているので、手を動かしながら解析を行っていくことができます。
まとめ
これは、完全に僕の本棚的なものなので、参考になるかどうかわかりませんが、参考になった方がいらっしゃれば幸いです。
pythonで一般化線形モデルと公式リファレンスの罠 - using statsmodels -
寒くなってきました。最近、pythonでデータの解析をすることにいそしんでおります。
Rでできることをpythonでやりたいなと思っていろいろ調べてみると、まぁなかなかできるようになっていなかったりするわけで、その辺を整備し始めたので、ここに書いていこうと思います。
一般化線形モデル（Generalized Linear Model）
一般化線形モデルは、回帰の問題ではよく利用される方法で、Rではglmという関数で使うことができます。有名なものはロジスティック回帰です。図的には、以下のような感じになります。
これをpythonで行うときには、多くの場合scikit-learnを使うような気がするのですが、scikit-learnのアウトプットは、統計的に結果を見たい場合に十分とはいえません。できればRのロジスティック回帰の結果みたいなの（下図）、、、を出力したいんです。
statsmodelsを使ってみよう。
そこで、そんな要望に答えるために、statsmodelsというモジュールが提供されています。どうもこれを使用すれば、Rのglm的なコトができるらしいと聞きつけて、やってみました。以下のコードをのせているので興味のある人は是非やってみて欲しいと思います。最後にRの対応するコードも付けておきました。使用しているデータはstatsmodelsの中のものですが、以下のリンクからもダウンロードできるようにしました。
https://dl.dropboxusercontent.com/u/2133906/sample.csv
※最後に注意事項も書いていますので、そこもできれば、読んでください。
## 一般化線形モデル
from __future__ import print_function
import numpy as np
import pandas as pd
from pandas import DataFrame, Series
from matplotlib import pyplot as plt
#glmを可能にするために必要なものなど
#statmodelsのGLMを使えるようにするモジュール
import statsmodels.api as sm
#GLMの中で用いる統計処理関数のインポート, chi**2の値などの算出に使用している
from scipy import stats
#Rのglmを使用可能にするための関数
import statsmodels.formula.api as smf
# ## OLS, GLM: Gaussian response data
#sampleデータのロード
#形式はpandas.DataFrame, naは削除する必要あり
df = sm.datasets.get_rdataset("Guerry", "HistData").data
df = df[['Lottery', 'Literacy', 'Wealth', 'Region']].dropna()
df.head()
#線形回帰モデルの構築, 最小自乗法を使う場合
#modelの構築
ols_model = 'Lottery ~ Literacy + Wealth + Region'
#モデルの構築
mod = smf.ols(formula=ols_model, data=df)
#fitに回帰した結果が入っているので、これをresに代入する
res = mod.fit()
#結果の参照
res.summary()
"""
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: Lottery No. Observations: 85
Model: GLM Df Residuals: 78
Model Family: Gaussian Df Model: 6
Link Function: identity Scale: 436.433771252
Method: IRLS Log-Likelihood: -375.30
Date: Fri, 07 Nov 2014 Deviance: 34042.
Time: 18:59:36 Pearson chi2: 3.40e+04
No. Iterations: 4
===============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
-------------------------------------------------------------------------------
Intercept 38.6517 9.456 4.087 0.000 20.118 57.186
Region[T.E] -15.4278 9.727 -1.586 0.113 -34.493 3.637
Region[T.N] -10.0170 9.260 -1.082 0.279 -28.167 8.133
Region[T.S] -4.5483 7.279 -0.625 0.532 -18.815 9.718
Region[T.W] -10.0913 7.196 -1.402 0.161 -24.195 4.013
Literacy -0.1858 0.210 -0.886 0.376 -0.597 0.225
Wealth 0.4515 0.103 4.390 0.000 0.250 0.653
===============================================================================
"""
#線形回帰モデルの構築, glmによるニュートン法を用いる場合
#modelの構築
glm_gauss_model = 'Lottery ~ Literacy + Wealth + Region'
#モデルの構築
mod = smf.glm(formula=glm_gauss_model, data=df, family=sm.families.Gaussian())
#fitに回帰した結果が入っているので、これをresに代入する
res = mod.fit()
#結果の参照
res.summary()
# ## GLM: Binomial response data
#
# ### Load data
#
# In this example, we use the Star98 dataset which was taken with permission
# from Jeff Gill (2000) Generalized linear models: A unified approach. Codebook
# information can be obtained by typing:
#一応データセットに関するディスクリプションを参照しておく。
print(sm.datasets.star98.NOTE)
"""
簡単に言えば、目的変数に「数学の成績が国の平均以上」だった人と、
「以下だった人」の数がたぶん学校毎に入っていて、それを学校に関係する
その他の変数で説明しましょうってことだと思います。
"""
#データセットの読み込み
star98 = sm.datasets.star98.load()
#pandas形式にしたい場合は
#star98 = sm.datasets.star98.load_pandas().data
"""
大事：
本来なら上記のロードだけで良いのだが、、、
pythonのstatsmodelsの通常のGLM = sm.GLMは
どういうわけか、intercept項を入れないまま回帰するので
sm.GLMを用いる場合には定数項を入れる必要があります。
"""
#data.exog = sm.add_constant(data.exog, prepend=False)
#datasetの形式を確認
type(star98.data)
# >>> numpy.core.records.recarray
#データセットをpandasのDataFrameにしておく。
star98 = DataFrame(star98.data)
#データセットの最初の5件を確認
star98.head()
#一般化線形モデルの構築
#2項回帰モデル、link:logit
"""
Binomialでモデリングする際には、
通常のRでは目的変数にcbind(成功数, 失敗数)を取るようにする。
これをpythonのstatsmodelsでは、~の前に以下のように記述する
'成功数 + 失敗数 ~ 説明変数'
"""
#モデルを定義する。
formula_binom = 'NABOVE + NBELOW ~ LOWINC+ PERASIAN + PERBLACK + PERHISP + PCTCHRT + PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF'
#成功割合の変数を作成する
star98['SUCCESS'] = star98['NABOVE'] / (star98['NABOVE'] + star98['NBELOW'])
star98['TOTAL'] = star98['NABOVE'] + star98['NBELOW']
# binomialの回帰モデルを構築する
glm_binom = smf.glm(formula=formula_binom, data=star98, family=sm.families.Binomial())
# res_binomにglmの結果を保存する
res_binom = glm_binom.fit()
# summaryを表示
res_binom.summary()
Generalized Linear Model Regression Results
================================================================================
Dep. Variable: ['NABOVE', 'NBELOW'] No. Observations: 303
Model: GLM Df Residuals: 282
Model Family: Binomial Df Model: 20
Link Function: logit Scale: 1.0
Method: IRLS Log-Likelihood: -2998.6
Date: Fri, 07 Nov 2014 Deviance: 4078.8
Time: 19:58:26 Pearson chi2: 4.05e+03
No. Iterations: 7
============================================================================================
coef std err z P>|z| [95.0% Conf. Int.]
--------------------------------------------------------------------------------------------
Intercept 2.9589 1.547 1.913 0.056 -0.073 5.990
LOWINC -0.0168 0.000 -38.749 0.000 -0.018 -0.016
PERASIAN 0.0099 0.001 16.505 0.000 0.009 0.011
PERBLACK -0.0187 0.001 -25.182 0.000 -0.020 -0.017
PERHISP -0.0142 0.000 -32.818 0.000 -0.015 -0.013
PCTCHRT 0.0049 0.001 3.921 0.000 0.002 0.007
PCTYRRND -0.0036 0.000 -15.878 0.000 -0.004 -0.003
PERMINTE 0.2545 0.030 8.498 0.000 0.196 0.313
AVYRSEXP 0.2407 0.057 4.212 0.000 0.129 0.353
PERMINTE:AVYRSEXP -0.0141 0.002 -7.391 0.000 -0.018 -0.010
AVSALK 0.0804 0.014 5.775 0.000 0.053 0.108
PERMINTE:AVSALK -0.0040 0.000 -8.450 0.000 -0.005 -0.003
AVYRSEXP:AVSALK -0.0039 0.001 -4.059 0.000 -0.006 -0.002
PERMINTE:AVYRSEXP:AVSALK 0.0002 2.99e-05 7.428 0.000 0.000 0.000
PERSPENK -1.9522 0.317 -6.162 0.000 -2.573 -1.331
PTRATIO -0.3341 0.061 -5.453 0.000 -0.454 -0.214
PERSPENK:PTRATIO 0.0917 0.015 6.321 0.000 0.063 0.120
PCTAF -0.1690 0.033 -5.169 0.000 -0.233 -0.105
PERSPENK:PCTAF 0.0490 0.007 6.574 0.000 0.034 0.064
PTRATIO:PCTAF 0.0080 0.001 5.362 0.000 0.005 0.011
PERSPENK:PTRATIO:PCTAF -0.0022 0.000 -6.445 0.000 -0.003 -0.002
============================================================================================
# ### 興味のある統計量を表示する方法
#試行回数
print("Total number of trials:" + str(star98[['NABOVE','NBELOW']].sum().sum(axis=1)))
#パラメータの表示
print('Parameters: ', res_binom.params)
#各パラメータに対するt値
print('T-values: ', res_binom.tvalues)
# ### Plots
# ロジスティック回帰の結果について図で確認する。
# プロットを通して、データの予測精度について分析する
# ・予測値と実測値のプロット
# ・予測値と残差のプロット
#観測値
y = star98['SUCCESS']
#推定値
yhat = Series(res_binom.mu)
#linear_predictor
y_lp = Series(np.log(res_binom.mu/(1-res_binom.mu)))
#残差
residuals = Series(res_binom.resid_pearson)
#dataframe:
dd = DataFrame({
'OBS':y,
'FITTED':yhat,
'LINEAR_PRED':y_lp,
'RESIDUALS':residuals
})
# 予測値と観測値の差を確認する
from ggplot import *
p_fitted_vs_obs = ggplot(aes(x='FITTED',y='OBS'),data=dd)+geom_point(colour="red")+geom_abline(intercept=0,slope=1)
# 予測値と残差の関係性を確認する
p_fitted_vs_residuals = ggplot(aes(x='FITTED',y='RESIDUALS'),data=dd)+geom_point(colour="red")+stat_smooth()
Rのコード
以下のようになります。よって、結果が一致しました！
> star98 = read.csv("sample.csv")
>
> glm_binom = glm(cbind(NABOVE,NBELOW) ~ LOWINC + PERASIAN + PERBLACK + PERHISP + PCTCHRT + PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF, data=star98, family="binomial")
> summary(glm_binom)
Call:
glm(formula = cbind(NABOVE, NBELOW) ~ LOWINC + PERASIAN + PERBLACK +
PERHISP + PCTCHRT + PCTYRRND + PERMINTE * AVYRSEXP * AVSALK +
PERSPENK * PTRATIO * PCTAF, family = "binomial", data = star98)
Deviance Residuals:
Min 1Q Median 3Q Max
-9.5217 -2.6505 -0.3357 2.1658 15.4200
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.959e+00 1.547e+00 1.913 0.0557 .
LOWINC -1.682e-02 4.340e-04 -38.749 < 2e-16 ***
PERASIAN 9.925e-03 6.014e-04 16.505 < 2e-16 ***
PERBLACK -1.872e-02 7.436e-04 -25.182 < 2e-16 ***
PERHISP -1.424e-02 4.339e-04 -32.818 < 2e-16 ***
PCTCHRT 4.917e-03 1.254e-03 3.921 8.81e-05 ***
PCTYRRND -3.580e-03 2.255e-04 -15.878 < 2e-16 ***
PERMINTE 2.545e-01 2.995e-02 8.498 < 2e-16 ***
AVYRSEXP 2.407e-01 5.714e-02 4.212 2.53e-05 ***
AVSALK 8.041e-02 1.392e-02 5.775 7.69e-09 ***
PERSPENK -1.952e+00 3.168e-01 -6.162 7.19e-10 ***
PTRATIO -3.341e-01 6.126e-02 -5.453 4.95e-08 ***
PCTAF -1.690e-01 3.270e-02 -5.169 2.36e-07 ***
PERMINTE:AVYRSEXP -1.408e-02 1.905e-03 -7.391 1.46e-13 ***
PERMINTE:AVSALK -4.005e-03 4.740e-04 -8.450 < 2e-16 ***
AVYRSEXP:AVSALK -3.906e-03 9.624e-04 -4.059 4.92e-05 ***
PERSPENK:PTRATIO 9.171e-02 1.451e-02 6.321 2.60e-10 ***
PERSPENK:PCTAF 4.899e-02 7.452e-03 6.574 4.89e-11 ***
PTRATIO:PCTAF 8.041e-03 1.499e-03 5.362 8.22e-08 ***
PERMINTE:AVYRSEXP:AVSALK 2.220e-04 2.989e-05 7.428 1.10e-13 ***
PERSPENK:PTRATIO:PCTAF -2.249e-03 3.490e-04 -6.445 1.15e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 34345.4 on 302 degrees of freedom
Residual deviance: 4078.8 on 282 degrees of freedom
AIC: 6039.2
Number of Fisher Scoring iterations: 4
公式リファレンスの罠。
githubにあるロジスティック回帰のコードが入っていますが、これはおそらく間違いです。
statsmodels/glm_formula.py at master · statsmodels/statsmodels · GitHub
どういうわけか、以下のコマンドでロジスティック回帰が実行できますということになっていますが...
## Generalized Linear Models (Formula)
# This notebook illustrates how you can use R-style formulas to fit Generalized Linear Models.
#
# To begin, we load the ``Star98`` dataset and we construct a formula and pre-process the data:
from __future__ import print_function
import statsmodels.api as sm
import statsmodels.formula.api as smf
star98 = sm.datasets.star98.load_pandas().data
formula = 'SUCCESS ~ LOWINC + PERASIAN + PERBLACK + PERHISP + PCTCHRT + PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF'
dta = star98[['NABOVE', 'NBELOW', 'LOWINC', 'PERASIAN', 'PERBLACK', 'PERHISP',
'PCTCHRT', 'PCTYRRND', 'PERMINTE', 'AVYRSEXP', 'AVSALK',
'PERSPENK', 'PTRATIO', 'PCTAF']]
endog = dta['NABOVE'] / (dta['NABOVE'] + dta.pop('NBELOW'))
del dta['NABOVE']
dta['SUCCESS'] = endog
# Then, we fit the GLM model:
mod1 = smf.glm(formula=formula, data=dta, family=sm.families.Binomial()).fit()
mod1.summary()問題があるのは、formulaの部分です。この形式では「目的変数」にSUCCESS（成功割合）のみを取ればよいということになっていますが、ロジスティック回帰の結果というのは「どれだけの回数試行を行ったか」に依存するため、成功の割合のみを回帰すれば良いというのは間違いです。実際、この結果を確認すると、以下のようになります。
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: SUCCESS No. Observations: 303
Model: GLM Df Residuals: 282
Model Family: Binomial Df Model: 20
Link Function: logit Scale: 1.0
Method: IRLS Log-Likelihood: -189.70
Date: Fri, 07 Nov 2014 Deviance: 380.66
Time: 20:11:04 Pearson chi2: 8.48
No. Iterations: 7
============================================================================================
coef std err z P>|z| [95.0% Conf. Int.]
--------------------------------------------------------------------------------------------
Intercept 0.4037 25.036 0.016 0.987 -48.665 49.472
LOWINC -0.0204 0.010 -1.982 0.048 -0.041 -0.000
PERASIAN 0.0159 0.017 0.910 0.363 -0.018 0.050
PERBLACK -0.0198 0.020 -1.004 0.316 -0.058 0.019
PERHISP -0.0096 0.010 -0.951 0.341 -0.029 0.010
PCTCHRT -0.0022 0.022 -0.103 0.918 -0.045 0.040
PCTYRRND -0.0022 0.006 -0.348 0.728 -0.014 0.010
PERMINTE 0.1068 0.787 0.136 0.892 -1.436 1.650
AVYRSEXP -0.0411 1.176 -0.035 0.972 -2.346 2.264
PERMINTE:AVYRSEXP -0.0031 0.054 -0.057 0.954 -0.108 0.102
AVSALK 0.0131 0.295 0.044 0.965 -0.566 0.592
PERMINTE:AVSALK -0.0019 0.013 -0.145 0.885 -0.028 0.024
AVYRSEXP:AVSALK 0.0008 0.020 0.038 0.970 -0.039 0.041
PERMINTE:AVYRSEXP:AVSALK 5.978e-05 0.001 0.068 0.946 -0.002 0.002
PERSPENK -0.3097 4.233 -0.073 0.942 -8.606 7.987
PTRATIO 0.0096 0.919 0.010 0.992 -1.792 1.811
PERSPENK:PTRATIO 0.0066 0.206 0.032 0.974 -0.397 0.410
PCTAF -0.0143 0.474 -0.030 0.976 -0.944 0.916
PERSPENK:PCTAF 0.0105 0.098 0.107 0.915 -0.182 0.203
PTRATIO:PCTAF -0.0001 0.022 -0.005 0.996 -0.044 0.044
PERSPENK:PTRATIO:PCTAF -0.0002 0.005 -0.051 0.959 -0.010 0.009
============================================================================================
回帰係数が異なっている分には実際、問題はなくて、問題があるのは「Yの推定値」に差があるかどうかです。それを確認するために、結果を確認します。まず、正しい結果はこのようになっています。
In [280]: yhat.head() Out[280]: 0 0.583312 1 0.751447 2 0.500583 3 0.685345 4 0.322510 dtype: float64
一方で間違っている結果は、このようになります。
In [279]: Series(mod1.mu).head() Out[279]: 0 0.578461 1 0.769802 2 0.449346 3 0.683891 4 0.247756 dtype: float64
先ほども説明した通り、公式リファレンスの罠ですね。下の結果は正しいロジスティック回帰の結果ではありません。この原因を作っているのは、先ほども言ったように成功割合のみでロジスティック回帰を実行したせいです。なぜなら、本来の成功回数Y(i)は、試行数N(i)、成功確率Pi(i)の2項分布に従います。しかしながら、Y(i)/N(i)の分布は当たり前ですが2項分布ではありません。確率変数変換をすれば自明です。なぜなら、P(i)の値域は0-1ですからね。結果的に、ロジスティック回帰は結果が違うことになってしまうというわけです。
注意）実際、成功割合のみを与えるようなロジスティック回帰が実行される場合には「全てのケースでの試行回数」が同じであるという仮定のもとで実行されていますので、結果は異なって当然です。
なので、pythonで行う際にはその点はご注意ください。
最後に
今回、pythonでglmを行う方法を考えましたが、まぁ実際Rでしようよっていうのが正直なところです。残念なことに笑
メモ@Python - sys.pathの問題「Extra/lib/python」という罠
夏休み以来、一度も投稿していませんでした。完全に、ブログ書くのも忘れて、研究にいそしんでおります．．．最近、Rからpythonに色々と移行しつつありまして、環境設定で、原因不明のエラーに苦しめられております。そろそろ、またいろんなlibraryを開拓したのでちょこちょこそういうのも載せていこうと思います。
python sys.path問題
ここ数日、pythonでnumpy, pandasなどをアップデートしたのに、アップデートしたものが、どういうわけだかimportすると昔のバージョンが使用される問題がありまして、かなり困っておりました。どうも、sys.pathに含まれている、以下のpathが悪さをしていることが分かりました。
sys.pathの中身
/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip /System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7 /System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin /System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac /System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages /System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python /System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk /System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old /System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload /Library/Python/2.7/site-packages
悪さしてるpath
'/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python'
なんだろこのpath...って思って、色々調べてみたのですが、こんな記事が出てきました...
Which Python · MacPython/wiki Wiki · GitHub
タイトルをみれば分かる通り...
--OS X の Pythonシステムを使うべきではない理由
ということです。
で、その中にこんなことが書かれていました。
OSX 10.9のシステムでは、numpy, scipy, matplotlibや、ウェブ系で便利なパッケージが入っています。そして、そのpathとして、"/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python"というものを採用しているようです。
その結果、新しいパッケージをインストールすると、
"/Library/Python/2.7/site-packages"
に本来ならインストールされて、pythonのpath内にこれは含まれているので、参照されるわけですが、このpathが"/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python"のpathよりも「下位」のものとしてリストに含まれています。なので、まずは/Extras/lib/pythonを先に参照してから、/site-packagesを読むため、2つのサイトに「numpy」が含まれていれば、先に読まれる方が優先されます。
そして、pip install numpyでインストールすると「site-packages」に入ってしまうため、いくらpipで最新版をインストールしても"/Extra"に入っているnumpy1.6.1が参照されます。などということで、標準搭載のpythonは非常に使いにくいということでした。原文は以下。
The interesting entry is "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python" because this entry is specific to system Python. Call this path "Extras/lib". On a standard OSX 10.9 machine this directory contains some Python packages including the numerical packages numpy, scipy and matplotlib and other useful web packages including twisted.
If you, the user, install new packages in system Python, the packages will get written to /Library/Python/2.7/site-packages. As you can see from the path listing above, your newly installed package will therefore be lower in the Python module path than the packages in Extras/lib. That means that if you try to install a new version of one of the packages in Extras/lib into system Python, system Python will ignore the version you installed and continue to use the version in Extras/lib.
Apple presumably chose this setup because they use Python for system tasks, and they want to be able to depend on working versions of the packages installed in Extras/lib. This is another reason to prefer another Python for your daily Python work.
相関の概念をさらにもう一歩深めるために。-ケンドールのタウを考えてみる-
昨日のエントリーで相関係数の有意性を確かめるためにどうすればいいのかということを書きました。
相関係数の有意性を確かめる方法について -相関係数について1歩踏み込む- - Data Science by R and Python
ですが、最後に書いたようにこの指標の有効性が検証できるのは「2つの変数が正規分布」に従うときに限られるということも最後に書きました。変数で正規分布に従わない変数はたくさんあります。身近な変数では「体重」は正規分布には従わないとされています。裾の重い分布、下の図のようになります。
URLBody Weight Distribution of Adult CT Scan Patients
この図、リンクを見ればわかる通りですが、右に引きづられています。この結果、分布は正規分布しないことがわかります。このようなデータ、変数同士の相関はどのようにして確認すればよいかという点を説明してみようと思います。
注意）以降の議論で使う「相関」というのは、これまでのような「直線的関係性」を意味していません。相関＝関連性を意味するものとします。
ケンドールのタウ（
）
ケンドールのという指標（ケンドールの順位相関係数）を取り上げます。この指標は変数X, Yの単調性の指標です。単調性というのは「直線的」な増加や減少を指していません。「Xが増加すれば、Yは増加する」「Yが減少すれば、Xも減少する」などの傾向を指します。それが「直線的な増加」や「直線的な減少」ではなくてもよいのです。指数的、対数的な関係であっても、単調でありさえすれば良いという便利な指標です。
ただし、どのような関係で変動するのかはわかりませんし、さらに「増加、減少」などが周期的に起こるような場合にも検出ができません。この点に付いては注意する必要があります。
しかしながら、この指標の優れているところは、相関係数とは異なり、2つの変数X,Yに特別な分布（正規分布など）を仮定しないという点です。このような手法を「ノンパラメトリック」と呼びます。
ケンドールのの詳しい説明
まず、ケンドールのの定義を書くと、
を同一の2変量分布に従う独立な組であるとします。このとき、ケンドールのタウは次のように定義されます。
さて、sgn{x}というのはxが正であれば1を返し、負であれば-1となるような関数です。この点を踏まえると、組の間にが成立すれば、増加傾向があり、逆に-1が成立するときには、減少の傾向があるということがわかります。これが定義です。
しかしながら、母集団のことは私たちにはわからないので、次のようなケンドールのタウの不偏推定量を用いてこの値をサンプルから推定するということを行います。なので、利用する際にはこちらを使います。
この値は実際、期待値を取ると上で定義したになることが示されます（簡単なので、興味のある人はやってみてください）。
このケンドールのタウという指標は、相関係数とにたような性質を持っています。例えば、変数XとYが独立であるとき、ケンドールの = 0が成立することが示されます。もちろん、逆は成立しません。相関が0 -> XとYは独立は成立しません。
ただし、ケンドールのタウが1となるような場合は相関係数の場合と状況が異なります。先ほども説明した通り、ケンドールのタウは「単調性」に関する指標なので、XとY=X**2のタウの値は1になります。一方で、相関係数は1にはなりません。
また、データのノイズが大きいと、ケンドールのタウは検出力が大幅に低下します。この点は注意が必要になります。
実際に実行する方法
これを実際に使うためのコマンドが以下になります。データについては、ニューヨークの気象データがRに初期装備されていたので、これを使います。まずは、データを図示するところから始めます。
head(airquality)
Ozone Solar.R Wind Temp Month Day
1 41 190 7.4 67 5 1
2 36 118 8.0 72 5 2
3 12 149 12.6 74 5 3
4 18 313 11.5 62 5 4
5 NA NA 14.3 56 5 5
#必要な変数は最初の4つ
df = airquality[,1:4]
#散布図とヒストグラムを書く
pairs3 <- function (data) {
oldpar <- par(no.readonly = TRUE) # 現在のグラフィックスパラメータ退避
on.exit(par(oldpar)) # (関数がエラー中断しても)パラメータ復帰
panel.hist <- function(x, ...) # ヒストグラムを描くための関数を定義
{ usr <- par("usr") # 現在のユーザー領域座標情報を得る
on.exit(par(usr)) # 関数終了時に usr パラメータ復帰
par(usr = c(usr[1:2], 0, 1.5) ) # ユーザー領域座標情報を変更
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)
}
pairs(data, panel=panel.smooth,
cex = 1.0, pch = 3, bg="light blue",
diag.panel=panel.hist, cex.labels = 2, font.labels=2)
}
#散布図 + HISTOGRAM
pairs3(df)
この図からわかるように、Ozoneと、Solarは正規分布に近くなさそうです。この2つを前回のエントリーの関数でqqplot, histogram重ね合わせた視覚的確認、さらにコルモゴロフ・スミルノフ検定を行ってみます。すると、以下の図と結果が得られます（結果はOzoneのみ確認します）。
> ks.test(na.omit(df[,1]),pnorm,mean=mean(na.omit(df[,1])),sd=sd(na.omit(df[,1]))) One-sample Kolmogorov-Smirnov test data: na.omit(df[, 1]) D = 0.148, p-value = 0.01243 alternative hypothesis: two-sided 警告メッセージ: In ks.test(na.omit(df[, 1]), pnorm, mean = mean(na.omit(df[, 1])), : コルモゴロフ・スミノフ検定において、タイは現れるべきではありません
ここで、警告メッセージが出されます。が、これは詳しくやり始めると大変なので無視しておきます。ks.testは値の中に「タイ（同じ数値）」があっては本当はいけないんです。
ただ、結果からもわかるように、5%有意で正規分布であるという帰無仮説を棄却することができ、かつ視覚的確認からも、qqnormの裾の方がかなり外れていることがわかります。よって、正規分布と仮定するのは適切ではないという結論を得ることができました。
変数が正規分布に従わないときの検定
このような場合、例えばOzoneとTempの関係性の数値を確認したいと思って「相関係数」を計算しても、この同時分布は前回も書きましたが「2変量正規分布」ではありませんので、検定をおこなうことはできません。なので、前回の手法で有意かどうかを確認することはできないということになります。
そこで、ケンドールのタウの登場です。ケンドールのタウを実際に計算すると以下のような値が得られます。
> cor(df$Ozone,df$Temp,use="pairwise",method="kendall") [1] 0.5862988
そして、この統計量が有意であるかどうかを調べるためには、前回同様検定を行う必要があります。このとき「検定」のためには以下の統計量を使います（この統計量は、漸近的に標準正規分布に従います。またKは上で求めたケンドールのタウです）。
この値を、実際に計算して検定するための関数が次の関数です。
> cor.test(df$Ozone,df$Temp,use="pairwise",method="kendall")
Kendall's rank correlation tau
data: df$Ozone and df$Temp
z = 9.1599, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
0.5862988 以上の結果より、「z = 9.1599, p-value < 2.2e-16」であるから、帰無仮説「ケンドールのタウは0」という仮説が成立する可能性は極めて低く棄却することができ、対立仮説「ケンドールのタウは0ではない」を採用することができます。よって、オゾン（Ozone）レベルと、気温（Temp）の間には「関連性」があることが認められました。
まとめ
さて、今回はケンドールのタウを説明してきました。おそらく、これまでの「相関係数」の考え方と随分離れており、取っ付きにくい方もいるかもしれません。ただ、正規分布しないデータに対して、正規分布だと決めつけて分析を行うのは、やはり危険ですし、少し難しくても別の手法を引き出してくる方が良いと思います。特に、ケンドールのタウはノンパラメトリックですから、ある意味どんな分布に従うデータに対しても使うことができます。
ただ、直感的に理解しづらい部分もありましたので、次回は「相関係数」とよく似ているけど、「ノンパラメトリック」な考え方である「スピアマンの順位相関係数」について書こうと思います。読んでくださってありがとうございました。
（告知）イベントを開催します。
統計の初心者の方（少しやったことある人も対象）向けに、データフェストというイベントを開催します。2日間でRの操作+実際の解析を行ってみようという試みです。興味のある方は、以下のリンクから詳細をご確認ください。（サイトが90年代後半感があるのは使用です笑）
http://datafest.jp/
p.s.
今回の、ブログごちゃごちゃしてる...
ブログで複雑な話するのはやっぱり無理があるのかな...
相関係数の有意性を確かめる方法について -相関係数について1歩踏み込む-
相関係数とは？
相関係数については、８月のエントリーを参照してもらうと良いかと思います。
「相関係数」ってなんですか？ -意味と利点と欠点をわかりやすく- - Data Science by R and Python
簡単に説明しておくと、相関というのは「2つの変数の間に存在する、直線的関係」を指しています。相関係数というのは、その直線度合いを0~1の数字で表したものと言えます。ただし、値には大小関係しかありません（相関係数0.8は、相関係数0.4のときよりも2倍相関が強いという意味ではないので注意してください）。
さて、データの分析・解析に関わらず、いくつかの変数があるデータでは、相関係数を計算して、変数の間に関係があるかをチェックするのが一般的です。前のエントリーにもあるようにこれは「直線的な関係」しか調べられないですし、色々と制約もあり、いちいち確認するのは面倒ですが、それでも使い勝手がいいのでみんな使っています。
さて、今日のモチベーションは相関係数を説明することではありません。「相関係数」を求めたときに、果たして求めた相関係数が「有意」であるか、「有意」ではないかに興味があるとします。つまり、求められた「相関係数」は統計的に意味がある数字であるのか、それとも「偶然」や「見かけ上」相関があるようにみえているだけで、本当は相関がない可能性もあるのか。この点を調べていきたいとします。つまり、「求めた相関係数が0である可能性があるか、ないか」を確認する方法についてです。
特に、相関係数を計算したあとに、その相関係数が「意味があるのか、ないのか」をチェックせずに進めてしまったりすると、問題になるかもしれませんし、「0.31の相関がありました！」のような、統計的に意味があることが保証されない報告をしてしまうことになるかもしれません。
なので、今回は求めた相関係数が有意かどうかを確かめる方法について考えていきます。
相関係数が統計的に有意（意味があるか）を検定する方法の前準備
まず、相関係数（ピアソンの積率相関係数）が有意かどうかを調べるためには、相関を計算する変数の同時分布が「2変量正規分布」に従うのかを調べておく必要があります（この仮定が成立しない場合は相関係数が有意かどうかを統計的に確かめることは難しいです）。
そのためには、次の事実を利用します。「2つの確率変数が、それぞれ正規分布（平均・分散は異なっても良い）に従うとき、その同時分布は正規分布に従う」。つまり、相関を確認したい2つの変数それぞれが正規分布に従うことを確認すればよいということになります。
直感的な方法（図で診断する！）
まず1つ目の簡単な方法は、「ヒストグラム」です。ヒストグラムの外形を見たときに、それが正規分布の外形と一致しているのかを確認します。そのためには、得られたデータから平均と分散を計算して、それらをパラメータにする理論的な正規分布をヒストグラムの上から書くのが良いでしょう。この方法はとても直感的で、おそらく確認の方法として粗くはありますが、強力です。
それから、もう1つはqqnormを用いる方法です。qqplotは、正規分布の累積分布関数を用います。細かい説明は、英語のwikipediaを確認してください。（Q–Q plot - Wikipedia, the free encyclopedia）。簡単に説明をしておくと、累積分布の25%点と75%点を直線で結び、その直線上にデータが乗るかを確認することで、正規分布に従っているかを確かめる方法です。
ここまで紹介してきた2つの方法を、手軽にできるようにするために、以下にRで作成した関数を添付しておきます。
x_norm = rnorm(1000,10,5)
x_unif = runif(1000)
range = function(val){
if((class(val) == 'numeric')| (class(val)=='vector')){
diff = (max(val) - min(val))/10
return(seq(min(val)-diff,max(val)+diff,(max(val)-min(val))/1000))
}else{
print("only one variable available")
}
}
norm_diag = function(x){
par(mfcol=c(1,2))
diff = (max(x) - min(x))/10
xrange = c(min(x)-diff,max(x)+diff)
hist(x,breaks=floor(length(x)/20),freq=F,xlim=xrange)
par(new=T)
points(range(x),dnorm(range(x),mean=mean(x), sd = sd(x)),type="l",xlim=xrange)
segments(x0=mean(x),y0=0,x1=mean(x),y1=dnorm(mean(x),mean(x),sd(x)),lty=2,col="red",lwd=1.5)
qqnorm(x,cex=0.5)
qqline(x,col="red",lty=2)
}
#変数が正規分布の場合
norm_diag(x_norm)
#結果はOKだと視覚的に確認できるはず
#変数が一様分布の場合
norm_diag(x_unif)
#結果はNOだと視覚的に確認できるはず
理詰めしたいときの検定
ですが、このままだと視覚的確認なので不安だという人のために、もう1つ統計的仮説検定を用いて正規分布に従うかどうかを検定する方法を紹介します。この検定は「コルモゴロフ・スミルノフ検定」と言います。仮説の設定をまずは確認しましょう。
ここで、注意が必要なのは「帰無仮説が棄却されれば、正規分布を仮定するのは適切ではない」という結論になることです。なので、正規分布を仮定しても、大きな問題にならないケースであることを示すには、P値が大きいことが必要です。
この検定の中身への言及はさけますが、簡単に言えば、経験分布関数と、理論的な分布関数の差の最大値を検定統計量として用いています。詳しくは、医学統計学の事典（丹後俊郎, 小西貞則 2010）を参考にしてください。この他、wikipediaなどにも記述があります。この検定をRで行うためには、以下のコマンドを実行すればOKです。簡単ですね。
#平均と分散をきちんと指定する必要があります #正規分布から発生させたデータの場合 > ks.test(x_norm,"pnorm",mean=mean(x_norm),sd=sd(x_norm)) # One-sample Kolmogorov-Smirnov test #data: x_norm #D = 0.02, <b><span style="color: #ff2600">p-value = 0.8185</span></b> #alternative hypothesis: two-sided #一様分布から発生させたデータの場合 > ks.test(x_unif,"pnorm",mean=mean(x_unif),sd=sd(x_unif)) One-sample Kolmogorov-Smirnov test data: x_unif D = 0.0599, p-value = 0.001522 alternative hypothesis: two-sided
結果を見れば、（上）正規分布から発生させた乱数では帰無仮説は棄却できませんが、（下）一様分布から発生したものは帰無仮説が棄却されてしまいます。よって、（上）は正規分布を仮定しても問題にはならなそうだということが推察され、（下）は正規分布を仮定するのは適切ではないという結論を得ます。
ちなみに、この他にも、検定手法はいろいろ提案されています。医学統計学の事典には、アンダーソン・ダーリング検定などが紹介されています。検定には、それぞれ特徴がありますので、他の検定手法では有意ではないなどの結果が出ることもあります。
※注意）ここで、注意ですが、「帰無仮説が棄却できなかった = 正規分布に従っている」は間違いです。帰無仮説が棄却できないというのは、帰無仮説が有意水準5%で棄却できないだけで、「変数が正規分布に従っている」と安易に結論づけてはいけません。
相関係数の検定
さて、ここまでの結果、相関を見たい変数が正規分布に従うと仮定しても問題ないとわかったとしましょう。すると、始めて相関係数が有意かどうかを検定することができます。
まずは、仮説の設定から始めます。
「帰無仮説：相関係数が0である」
「対立仮説：相関係数は0ではない」
ここで、「相関係数が0である」というのは、確率変数には正規分布を仮定していますから、独立性の検定と同等です。これを検定するためには、2つの方針があります。1つ目は尤度比検定です。もう1つが、帰無仮説の元での標本相関係数の分布を求めるというやり方です。今回は後者を考えます。
さて、帰無仮説の元での標本相関係数の分布は、t分布に従います。この証明は、非常に難解です。興味のある人は以下の資料をご覧下さい。おそらく、ほとんどの人がやらないでもいいんじゃないでしょうか．．．
Pearson product-moment correlation coefficient - Wikipedia, the free encyclopedia
Correlation Coefficient--Bivariate Normal Distribution -- from Wolfram MathWorld
こんな難しいことを知らなくても、Rを使えば簡単に検定を行うことができます。2つの正規分布に従う確率変数の間に相関があるかないかは、次のコマンドで実行できます。
#無相関 > cor.test(x_norm, rnorm(1000)) #Pearson's product-moment correlation # #data: x_norm and rnorm(1000) #t = 0.097, df = 998, p-value = 0.9227 #alternative hypothesis: true correlation is not equal to 0 <b><span style="color: #ff2600">#95 percent confidence interval: # -0.05893338 0.06505162</span></b> #sample estimates: # cor #0.003070921 #相関あり > cor.test(x_norm, x_norm + rnorm(1000,0,8)) #Pearson's product-moment correlation # #data: x_norm and x_norm + rnorm(1000, 0, 8) #t = 20.5626, df = 998, p-value < 2.2e-16 #alternative hypothesis: true correlation is not equal to 0 <span style="color: #ff2600"><b>#95 percent confidence interval: # 0.5004476 0.5876367</b></span> #sample estimates: # cor #0.5455164
これを確認して、有意であれば95%信頼区間が0をまたがない出力を得ます。前者は0を跨いでしまっているので、有意ではないですね。一方で、後者は0をまたいでないので、有意に相関があることがわかります。
まとめ
今回は、相関係数が有意であるかどうかの検定について書きました。相関係数が有意であるかどうかを統計的に確認するのは、様々な仮定が必要です。そして、満たさないデータもたくさんあります。ですが、満たすデータもあるはずですし、そのような場合には、上記のような検定は効果を発揮すると思います。
でも、満たさない場合も扱いたい！そんな人もいらっしゃるかもしれません。次回は、正規分布に従わない場合に変数同士の関係性（相関ではなく）を確認する方法の1つである「スピアマンの順位相関」について、それから見かけ上の相関（疑似相関）を取り除いて純粋に2つの変数の関係性を確認するための考え方である「偏相関」について書きたいと思います。
（告知）イベントを開催します。
統計の初心者の方（少しやったことある人も対象）向けに、データフェストというイベントを開催します。2日間でRの操作+実際の解析を行ってみようという試みです。興味のある方は、以下のリンクから詳細をご確認ください。（サイトが90年代後半感があるのは使用です笑）
http://datafest.jp/
PypeRでpythonからRを使う -- データ解析をpythonで！
PythonからRを使うための方法．
pythonからRを使うためには、pythonにpypeRというパッケージ？をインストールする必要があります。こちらです：PypeR
$pip install PypeR
で導入して使ってみてください。
実際にPypeRを使ってみます。
PypeRを使ってみた結果は以下から参照できます。ipython notebookが便利なので、これからはここにプログラムを書くのをやめて、こんな感じで公開するのもいいかな。。。と思ったりしています。
いや、本当に便利なんですnotebookが。そういえば、Rpubsとか、どうなんでしょう？使いやすいんでしょうか？PypeRを使ってみた結果を解説付きでこちらで公開しました。
http://nbviewer.ipython.org/github/tomoshige/python_blog/blob/master/pypeR/pypeR_sample.ipynb
まとめ
うーん、一応pythonでRは使えるけど、やっぱりまだまだ不便な部分もあるなぁと。。。慣れれば楽になるんですかね？特にいちいちRに渡して、結果を戻したりするのが一手間かかかるんですよね。やはり、pythonはデータ取得・データ加工に特化させて、Rで解析・解釈するのがいいんでしょうかね。でも、全部pythonで完結できるようにした方が後々いいんじゃないかなとか思ったりもしますね。
補足
さて、そろそろ、本格的にscikit-learn使い始めようとおもっているのと、最近クローリングとスクレイピングの勉強をしてるのですが、なかなか進みません（泣。。。おなかすいた...
pythonでggplotを使おう！ -pythonによるデータ分析（２）-
ggplot2をpythonで使う。
Rを使っている人にとっては、ggplot2はおなじみな感じがしますが、pythonでもこれ使えるといいなということで、pythonで使える方法を探しておりましたら、"ggplot"はpythonでもあるんですね。
ggplot 0.6.5 : Python Package Index
ということで、これを積極的に使っていくためにセットアップします。あんまり日本語の記事もないですし。ちなみに、Rではggplot2でしたが、pythonではggplotです。
依存パッケージなど
ggplotは以下のパッケージに依存しています、pythonでのデータ解析で、使うモノばかりなのでインストールしている人も多いかと思いますが、ない人はインストールをしてください。
- matplotlib
- pandas
- numpy
- scipy
- statsmodels
- patsy
インストール
インストールは、いつも通りpipを使ってできます。
$ pip install ggplot
基本操作
基本操作はRとほとんど変わらないので、気にすることはないと思います。これで、pythonでggplotができるようになりました。
