# この企画は、雑誌や教科書にでているグラフをSageで再現し、 # グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。 #
## 前回に続き、データ解析のための統計モデリング入門 # （以下、久保本と書きます） # の第6章の例題をSageを使って再現してみます。 #
## 数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるのが大きな特徴です。 # この機会にSageを分析に活用してみてはいかがでしょう。 #
# # # ## 最初に必要なライブラリーやパッケージをロードしておきます。 #
# # In[1]: # RとPandasのデータフレームを相互に変換する関数を読み込む # Rの必要なライブラリ r('library(ggplot2)') r('library(jsonlite)') # python用のパッケージ import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns # statsmodelsを使ってglmを計算します import statsmodels.formula.api as smf import statsmodels.api as sm from scipy.stats.stats import pearsonr get_ipython().run_line_magic('matplotlib', 'inline') # jupyter用のdisplayメソッド from IPython.display import display, Latex, HTML, Math, JSON # sageユーティリティ load('script/sage_util.py') # Rユーティリティ load('script/RUtil.py') # ## 6章の最初に二項分布を使った回帰分析の例がでてきます。ポアソン分布では上限のないカウントデータでしたが、 # 上限のある場合には、二項分布を使うのだと説明がありました。 # （自然界ではこのように上限のあるカウントデータが多いので、とても参考になりました） #
# # # ## 二項分布用のデータは、サポートページにあるdata4a.csvです。 # これを読み込んで、データの性質と分布を表示します。 #
## 直感的に施肥を施す(f=T)ことでサイズxが大きくなっているのが分かります。 # また、サイズと種子数の右上がりの関係が見られます。（種子数の上限は8です。） #
# # In[2]: # 6章のデータを読み込む d = pd.read_csv('data/data4a.csv') d.head() # In[3]: d.describe() # In[4]: # F別の分布をみる sns.lmplot('x', 'y', data=d, hue='f', fit_reg=False ) plt.show() # ## ここで、Sageのプロット関数を使って二項分布の確率分布がどのような形になっているか # 久保本の図6.3と同じ条件でブロットしてみます。 #
## 二項分布を # $$ # p(y | N, q) = \binom{N}{y} q^y (1 -q)^{N-y} # $$ # 以下の様に_p関数で定義します。 #
## Sageでプロットしてみます。Sageではグラフ（Graphics）にグラフを足すことで、重ね合わせのプロットができます。 # とても直感的にグラフを書くことができます。（ggplotも同じです） #
# # In[5]: # 二項分布を定義 def _p(q, y, N): return binomial(N, y)*q^y*(1-q)^(N-y) # In[6]: g = Graphics() g += list_plot([_p(0.8, y, 8) for y in (0..8)], plotjoined=True, rgbcolor="red", legend_label ="q=0.8") g += list_plot([_p(0.3, y, 8) for y in (0..8)], plotjoined=True, rgbcolor="green", legend_label ="q=0.3") g += list_plot([_p(0.1, y, 8) for y in (0..8)], plotjoined=True, rgbcolor="blue", legend_label ="q=0.1") g.show(figsize=4) # ## ロジスティック関数 # $$ # q_i = logistic(z_i) = \frac{1}{1 + exp(-z_i)} # $$ # の関数を_logisticとして定義し、その分布をSageを使ってプロットします。 #
## ロジスティック関数のような曲線を持つデータの場合には、二項分布の確率分布があると予測して回帰分析を行います。 # リンク関数は、ロジスティック関数の逆関数であるロジットリンク関数を指定します。 #
## ロジット関数は、以下の様に表されます(式a)。 # $$ # logit(q_i) = log \frac{q_i}{1 - q_i} # $$ #
# # In[7]: # ロジスティック関数の定義 def _logistic(z): return 1/(1 + exp(-z)) # In[8]: # ロジスティック曲線 plot(_logistic(x), [x, -6, 6], figsize=4, legend_label='$q=\\frac{1}{1+exp(-z)}$') # ## 残念ながら、statmodelsではcbindに相当する二項分布の解析方法が分からず、 # Rを使って計算することにしました。 #
## Sageの中からRの関数を使うことができるので、計算を中断することなく進めることができます。 #
## 回帰の結果、$\beta_1 = -19.536, \beta_2 = 1.95, \beta_3=2.02$と求まりました。 #
# # In[9]: # statmodelsでcbind(y, N -y)の部分を表現する方法が見つからなかった # 0, 1の場合には、statmodelsでも解析可能です。 # 残念ですが、Rを使って処理します。 PandaDf2RDf(d, "d") # In[10]: r('glm(cbind(y, N - y) ~ x + f, data=d, family=binomial)') # ### statmodelsで二項分布を解析 # 以下のサイトにstatmodelsで二項分布を解析する方法が紹介されていました。 # - [Data Science by R and Python](http://tomoshige-n.hatenablog.com/entry/2014/11/07/202201) # # cbind(A, B)の部分は、A + Bと表現すればよいことが分かりました。 # In[11]: # statmodelsで解析 # N-yのカラムを追加 N = 8 d['NmY'] = N - d.y fit = smf.glm('y + NmY ~ x + f', data=d, family=sm.families.Binomial()).fit() fit.summary2() # ## Sageを使って解析結果、施肥あり$logistic(-19.536+1.952x + 2.022)$、施肥なし$logistic(-19.536+1.952x)$をプロットしました。 #
# # In[12]: dT = d[d.f == 'T'] dC = d[d.f == 'C'] dT_pts = list_plot(zip(dT.x, dT.y), zorder=2) dC_pts = list_plot(zip(dC.x, dC.y), zorder=2, rgbcolor='red') dT_plt = plot(_logistic(-19.536+1.952*x + 2.022)*8, [x, 7, 13]) dC_plt = plot(_logistic(-19.536+1.952*x)*8, [x, 7, 13], rgbcolor='red') (dT_pts+dT_plt+dC_pts+dC_plt).show(figsize=4) # ## 久保本のすごいところは、単に回帰をすることで終わらず、その解釈を説明しているところです。 # リンク関数がロジット関数である(式a)であることから # $$ # log \frac{q_i}{1 - q_i} = \beta_1 + \beta_2 x_i + \beta_3 f_i # $$ # となり、両辺の指数を取ると、 # $$ # \frac{q_i}{1 - q_i} = exp(\beta_1 + \beta_2 x_i + \beta_3 f_i) # = exp(\beta_1) exp(\beta_2 x_i) exp(\beta_3 f_i) # $$ # となります。$\frac{q_i}{1 - q_i}$がオッズと呼ばれる量です。 #
## 施肥の有無による違いは、exp(2.02)=7.5倍であると推定されました。 #
# # ## 久保本の6章で感動したのは、割り算を使わない方法です。 # 私のようなものは、割り算を使って正規化してしまいますが、 # このように割り算をしないで、回帰分析する方法としてオフセット項の使い方を # 紹介しています。 #
## オフセット用のデータには、data4b.csvを使います。 # このデータは、 #
# オフセット項は、offset=np.log(オフセットのカラム)のように指定します。ここでは、オフセット値Aのlogを # 取った値が使われていますが、これは以下のような理由からです。 #
## 人口密度を以下の様に表し、 # $$ # \frac{平均個体数\lambda_i}{A_i} = 人口密度 # $$ # 平均個体数と明るさに指数関数的な関係があると仮定すると、 # $$ # \lambda_i = A_i \times 人口密度 = A_i exp(\beta_1 + \beta_2 x_i) # $$ # ここで、面積を掛けている部分にlogを使うことで # $$ # \lambda_i = exp(\beta_1 + \beta_2 x_i + log A_i) # $$ # と変形することができます。つまり、$log A_i$だけずれた値（オフセット項）、久保本では「げた」をはかせたと考えると # よいとありました。 #
# # In[13]: # オフセット用のデータを読み込む d2 = pd.read_csv('data/data4b.csv') d2.head() # In[14]: # Aの値で色を変えてプロット sns.lmplot('A', 'y', data=d2, hue='x', fit_reg=False, legend=False) plt.show() # ## 推定された値をプロットすると以下のようになります。 # データよりもはっきりと明るさによる影響がきれいに分かります。（そのようにモデルを作ったので当たり前なのですが） #
# # In[15]: fit = smf.glm('y ~ x', data=d2, offset=np.log(d2.A), family=sm.families.Poisson(link=sm.families.links.log)).fit() fit.summary2() # In[16]: d2['pred'] = fit.predict() # xの値で色を変えてプロット sns.lmplot('A', 'pred', data=d2, hue='x', fit_reg=False, legend=False) plt.show() # ## Sageにも確率分布を扱う関数RealDistirubtionがVer.5から導入されました。（まだ一部の関数のサポートですが） # これを使って、正規分布の確率密度関数をプロットしていました。 #
# # In[17]: # Sage ver.5から導入されたRealDistributionを使う T1 = RealDistribution('gaussian', 1, seed=101) T2 = RealDistribution('gaussian', 3, seed=101) # In[18]: T1_plt = plot(lambda x : T1.distribution_function(x), [x, -4, 4], rgbcolor="red") T2_plt = plot(lambda x : T2.distribution_function(x), [x, -4, 4], rgbcolor="green") T3_plt = plot(lambda x : T1.distribution_function(x-2), [x, -4, 4], rgbcolor="blue") (T1_plt + T2_plt + T3_plt).show(figsize=4) # ## 確率密度関数から指定された範囲に含まれる確率は、その範囲の面積になります。（その範囲で積分した値） # cum_distribution_functionは、マイナス無限大から指定された値までの確立を計算するので、 # この関数の差を計算することで、指定範囲の確率（面積）を求めることができます。 #
# # In[19]: # xが1.2から1.8となる確率は、cum_distribution_functionの差で計算で求まる T1.cum_distribution_function(1.8) - T1.cum_distribution_function(1.2) # ## 久保本に従って正規分布の尤度をと最尤尤度を求めてみます。 # $y_iがy_i - 0.5\Delta y \leq y \leq y_i + 0.5\Delta y$である確率は、 # $$ # \begin{eqnarray} # L(\mu, \sigma) & = & \prod_i p(y_i | \mu, \sigma) \Delta y \\ # & = & \prod_i \frac{1}{\sqrt{2 \pi \sigma^2}} exp\left\{ \frac{(y_i - \mu)^2}{2\sigma^2} \right\} \Delta y \\ # & = & \prod_i (2 \pi \sigma^2)^{-\frac{1}{2}} exp\left\{ \frac{(y_i - \mu)^2}{2\sigma^2} \right\} \Delta y # \end{eqnarray} # $$ # 両辺をlogを取って # $$ # \begin{eqnarray} # log L(\mu, \sigma) & = & -\frac{1}{2} \sum_i log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - \mu)^2 + \sum_i log(\Delta y) \\ # & = & -\frac{1}{2} N log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - \mu)^2 + N log(\Delta y) # \end{eqnarray} # $$ # $\Delta y$の項を無視すると、 # $$ # log L(\mu, \sigma) = -\frac{1}{2} N log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - \mu)^2 # $$ # となります。これを$\mu$で偏微分し、最小となるのは、$y = \mu$と求まります。（最小自乗法の解と同じ） #
# # ## どんなときにガンマ分布を使うのかGoogleでみると、正の値で、寿命や待ち時間など経過時間に無関係一定値の場合に、 # 使われるとありました。 #
## ガンマ分布の確立密度関数は、以下の様に定義されます。 # $$ # p(y\,|\,s, r) = \frac{r^s}{\Gamma(s)}y^{s-1}\, exp(-ry) # $$ # ここでsはshapeパラメータ、rはrateパラメータと呼ばれています。 # # ガンマ分布の平均は$s/r$、分散は$s/r^2$となります。 #
# # In[20]: # ガンマ分布の確率密度分布を定義 def _p(y, s, r): return r^s/gamma(s)*y^(s-1)*e^(-r*y) # In[21]: y = var('y') plt1 = plot(_p(y,1,1), [y,0,5], rgbcolor="red", legend_label ="r=s=1") plt2 = plot(_p(y,5,5), [y,0,5], rgbcolor="green", legend_label ="r=s=5") plt3 = plot(_p(y,0.1,0.1), [y,0,5], rgbcolor="blue", legend_label ="r=s=0.1") (plt1+plt2+plt3).show(figsize=4, ymax=1) # ## 例題では、ある個体の花の重量を$y_i$が平均$\mu_i$のガンマ分布にしたがっていると仮定します。 # $$ # \mu_i = A x_i^b # $$ # で、A=exp(a)とすると、 # $$ # \mu_i = exp(a)x_i^b = exp(a + log x_i^b) = exp(a + b log x_i) # $$ # と表され、リンク関数は、expの逆関数のlogとなり、 # $$ # log \mu_i = a + b log x_i # $$ # となります。 #
## ガンマ分布の例題データは、RData形式なので、Rで解析します。family=Gamma(link="log")でリンク関数logのガンマ分布を指定します。 #
## glmを使うと複雑な分布の回帰が簡単に計算できます。しかし、どのモデルが良いかは、 # プロット結果を見ながら判断するのがよいと実感しました。人間の目で分からない違いは、尤度やAICを使うのが現実できなのでは？ #
# # In[22]: # ガンマ分布用のデータを使って回帰分析 r('load("data/d.RData")') r('glm(y ~ log(x), family=Gamma(link="log"), data=d)') # In[23]: # 結果の図化、stat_smoothのglmを使う場合、method.argsにfamilyをセットするのに変わったみたい。 graph = preGraph("images/fig-6.13.pdf") r('p <- ggplot(data=d, aes(x=x, y=y)) + geom_point() + stat_smooth(method=glm, formula=y ~ log(x), method.args = list(family=Gamma(link="log")))') r('plot(p)') postGraph(graph) # ## おまけ # Rの代わりにPythonでガンマ分布のGLM例題を試してみます。 # In[24]: # Rからsageにデータを取り込みます。 d = RDf2PandaDf("d") # 結果のプロットを簡単にするために、xでソートしたデータをdにセット d = d.sort_values(by=['x']) # 分布をプロット d.plot(kind='scatter', x='x', y='y') plt.show() # In[25]: fit = smf.glm('y ~ np.log(x)', data=d, family=sm.families.Gamma(link=sm.families.links.log)).fit() fit.summary2() # In[26]: # 推定結果をd.predに追加して、プロット d['pred'] = fit.predict() # 分布をプロット ax = d.plot(kind='scatter', x='x', y='y') d.plot(kind='line', x='x', y='pred', ax=ax) plt.show() # In[ ]: