sage/データ解析のための統計モデリング入門勉強ノート第6章
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2014/04/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/43/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう：データ解析のための統計モ...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon....
（以下、久保本と書きます）
の第6章の例題をSageを使って再現してみます。
数式処理システムSageのノートブックは、計算結果を表示する...
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#v6dc0e11]
最初に必要なライブラリーやパッケージをロードしておきます。
sageへの入力:
#pre{{
# Rの必要なライブラリ
#r("install.packages('ggplot2')")
r('library(ggplot2)')
r('library(jsonlite)')
# python用のパッケージ
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from ggplot import *
# statsmodelsを使ってglmを計算します
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
** 上限のあるカウントデータの回帰分析 [#e8efce13]
6章の最初に二項分布を使った回帰分析の例がでてきます。ポア...
上限のある場合には、二項分布を使うのだと説明がありました。
（自然界ではこのように上限のあるカウントデータが多いので...
*** 二項分布用のデータについて [#o136f767]
二項分布用のデータは、サポートページにあるdata4a.csvです。
これを読み込んで、データの性質と分布を表示します。
直感的に施肥を施す(f=T)ことでサイズxが大きくなっているの...
また、サイズと種子数の右上がりの関係が見られます。（種子...
sageへの入力:
#pre{{
# 6章のデータを読み込む
d = pd.read_csv(DATA+'data4a.csv')
d.head()
}}
#pre{{
N y x f
0 8 1 9.76 C
1 8 6 10.48 C
2 8 5 10.83 C
3 8 6 10.94 C
4 8 1 9.37 C
[5 rows x 4 columns]
}}
sageへの入力:
#pre{{
d.describe()
}}
#pre{{
N y x
count 100 100.000000 100.000000
mean 8 5.080000 9.967200
std 0 2.743882 1.088954
min 8 0.000000 7.660000
25% 8 3.000000 9.337500
50% 8 6.000000 9.965000
75% 8 8.000000 10.770000
max 8 8.000000 12.440000
[8 rows x 3 columns]
}}
sageへの入力:
#pre{{
ggplot(d, aes(x='x', y='y', color='f')) + geom_point()
}}
#pre{{
<ggplot: (31457393)>
}}
sageへの入力:
#pre{{
ggsave('fig-6.2.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-6.2.png);
*** 二項分布 [#cfca1248]
ここで、Sageのプロット関数を使って二項分布の確率分布がど...
久保本の図6.3と同じ条件でブロットしてみます。
二項分布を
$$
p(y | N, q) = \begin{pmatrix} N \\ y \end{pmatrix} q^y (1...
$$
以下の様に_p関数で定義します。
Sageでプロットしてみます。Sageではグラフ（Graphics）にグ...
とても直感的にグラフを書くことができます。（ggplotも同じ...
sageへの入力:
#pre{{
# 二項分布を定義
def _p(q, y, N):
return binomial(N, y)*q^y*(1-q)^(N-y)
}}
sageへの入力:
#pre{{
g = Graphics()
g += list_plot([_p(0.8, y, 8) for y in (0..8)], plotjoine...
g += list_plot([_p(0.3, y, 8) for y in (0..8)], plotjoin...
g += list_plot([_p(0.1, y, 8) for y in (0..8)], plotjoin...
g.show(figsize=5)
}}
&ref(sage0.png);
*** ロジスティック関数 [#t5a0280e]
ロジスティック関数
$$
q_i = logistic(z_i) = \frac{1}{1 + exp(-z_i)}
$$
の関数を_logisticとして定義し、その分布をSageを使ってプロ...
ロジスティック関数のような曲線を持つデータの場合には、二...
リンク関数は、ロジスティック関数の逆関数であるロジットリ...
ロジット関数は、以下の様に表されます(式a)。
$$
logit(q_i) = log \frac{q_i}{1 - q_i}
$$
sageへの入力:
#pre{{
# ロジスティック関数の定義
def _logistic(z):
return 1/(1 + exp(-z))
}}
sageへの入力:
#pre{{
# ロジスティック曲線
plot(_logistic(x), [x, -6, 6], figsize=5, legend_label='$...
}}
&ref(sage1.png);
*** パラメータ推定 [#c2e896dd]
残念ながら、statmodelsではcbindに相当する二項分布の解析方...
Rを使って計算することにしました。
Sageの中からRの関数を使うことができるので、計算を中断する...
回帰の結果、\(\beta_1 = -19.536, \beta_2 = 1.95, \beta_3=...
sageへの入力:
#pre{{
# statmodelsでcbind(y, N -y)の部分を表現する方法が見つか...
# 0, 1の場合には、statmodelsでも解析可能です。
# 残念ですが、Rを使って処理します。
PandaDf2RDf(d, "d")
}}
sageへの入力:
#pre{{
r('glm(cbind(y, N - y) ~ x + f, data=d, family=binomial)')
}}
#pre{{
Call: glm(formula = cbind(y, N - y) ~ x + f, family = bi...
Coefficients:
(Intercept) x fT
-19.536 1.952 2.022
Degrees of Freedom: 99 Total (i.e. Null); 97 Residual
Null Deviance: 499.2
Residual Deviance: 123 AIC: 272.2
}}
*** 結果の図化 [#z6399d53]
ggplot2のstat_smoothで結果の曲線を表示しようとしたのです...
sageへの入力:
#pre{{
# glmの結果を使って曲線を表示できませんでした！stat_smoot...
graph = preGraph("fig-6.7.png")
#r('p <- ggplot(d, aes(x=x, y=y)) + geom_point() + facet_...
r('p <- ggplot(d, aes(x=x, y=y)) + geom_point() + facet_w...
r('plot(p)')
postGraph(graph)
}}
&ref(fig-6.7.png);
そこで、Sageのプロット機能を使って結果を表示してみました。
sageへの入力:
#pre{{
# 代わりにSageでプロットしてみました。
dT = d[d.f == 'T']
dC = d[d.f == 'C']
dT_pts = list_plot(zip(dT.x, dT.y))
dC_pts = list_plot(zip(dC.x, dC.y), rgbcolor='red')
dT_plt = plot(_logistic(-19.536+1.952*x + 2.022)*8, [x, 7...
dC_plt = plot(_logistic(-19.536+1.952*x)*8, [x, 7, 13], r...
(dT_pts+dT_plt+dC_pts+dC_plt).show(figsize=5)
}}
&ref(sage2.png);
*** 結果の解釈 [#g22c51be]
久保本のすごいところは、単に回帰をすることで終わらず、そ...
リンク関数がロジット関数である(式a)であることから
$$
log \frac{q_i}{1 - q_i} = \beta_1 + \beta_2 x_i + \beta_3...
$$
となり、両辺の指数を取ると、
$$
\frac{q_i}{1 - q_i} = exp(\beta_1 + \beta_2 x_i + \beta_3...
= exp(\beta_1) exp(\beta_2 x_i) exp(\beta_3 f_i)
$$
となります。\(\frac{q_i}{1 - q_i}\)がオッズと呼ばれる量で...
施肥の有無による違いは、exp(2.02)=7.5倍であると推定されま...
** 割り算を使わない方法はすごい！ [#l375b87e]
久保本の6章で感動したのは、割り算を使わない方法です。
私のようなものは、割り算を使って正規化してしまいますが、
このように割り算をしないで、回帰分析する方法としてオフセ...
紹介しています。
オフセット用のデータには、data4b.csvを使います。
このデータは、
- 調査地iの面積を\(A_i\)
- 調査地iの明るさを\(x_i\)
- 調査地iにおける植物個体数\(y_i\)
の記録です。横軸にA, 縦軸にyを取りデータをプロットすると...
面積が大きいほど植物個体数が多く、xが大きい（明るい）程植...
見て取れます。
オフセット項は、offset=np.log(オフセットのカラム)のように...
取った値が使われていますが、これは以下のような理由からで...
人口密度を以下の様に表し、
$$
\frac{平均個体数\lambda_i}{A_i} = 人口密度
$$
平均個体数と明るさに指数関数的な関係があると仮定すると、
$$
\lambda_i = A_i \times 人口密度 = A_i exp(\beta_1 + \beta...
$$
ここで、面積を掛けている部分をlogを使うことで
$$
\lambda_i = exp(\beta_1 + \beta_2 x_i + log A_i)
$$
と変形することができます。つまり、\(log A_i\)だけずれた値...
よいとありました。
sageへの入力:
#pre{{
# オフセット用のデータを読み込む
d2 = pd.read_csv(DATA+'data4b.csv')
d2.head()
}}
#pre{{
y x A
0 57 0.68 10.3
1 64 0.27 15.6
2 49 0.46 10.0
3 64 0.45 14.9
4 82 0.74 14.0
[5 rows x 3 columns]
}}
sageへの入力:
#pre{{
ggplot(d2, aes(x='A', y='y', color='x')) + geom_point()
}}
#pre{{
<ggplot: (40327389)>
}}
sageへの入力:
#pre{{
ggsave('fig-6.10-A.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-6.10-A.png);
*** 推定値 [#b4892baf]
推定された値をプロットすると以下のようになります。
データよりもはっきりと明るさによる影響がきれいに分かりま...
sageへの入力:
#pre{{
fit = smf.glm('y ~ x', data=d2, offset=np.log(d.A), famil...
fit.summary2()
d2['pred'] = fit.predict()
}}
sageへの入力:
#pre{{
# 予測値のプロット
ggplot(d2, aes(x='A', y='pred', color='x')) + geom_point()
}}
#pre{{
<ggplot: (39619317)>
}}
sageへの入力:
#pre{{
ggsave('fig-6.10-C.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-6.10-C.png);
*** オフセット項を使わないと [#k92b37cd]
オフセット項を使わないと面積による影響が回帰分析に含まれ...
sageへの入力:
#pre{{
# 人口密度Aをオフセットを使わない場合
fit_2 = smf.glm('y ~ x', data=d2, family=sm.families.Pois...
fit_2.summary2()
d2['pred2'] = fit_2.predict()
}}
sageへの入力:
#pre{{
ggplot(d2, aes(x='A', y='pred2', color='x')) + geom_point()
}}
#pre{{
<ggplot: (40798181)>
}}
sageへの入力:
#pre{{
ggsave('fig-6.10-D.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-6.10-D.png);
** 正規分布と尤度 [#m8c6dfd3]
Sageにも確率分布を扱う関数RealDistirubtionがVer.5から導入...
これを使って、正規分布の確率密度関数をプロットしていまし...
sageへの入力:
#pre{{
# Sage ver.5から導入されたRealDistributionを使う
T1 = RealDistribution('gaussian', 1, seed=100)
T2 = RealDistribution('gaussian', 3, seed=100)
}}
sageへの入力:
#pre{{
T1_plt = plot(lambda x : T1.distribution_function(x), [x,...
T2_plt = plot(lambda x : T2.distribution_function(x), [x,...
T3_plt = plot(lambda x : T1.distribution_function(x-2), [...
(T1_plt + T2_plt + T3_plt).show(figsize=5)
}}
&ref(sage3.png);
*** 確率密度関数から確率を得る [#u6ff9651]
確率密度関数から指定された範囲に含まれる確率は、その範囲...
cum_distribution_functionは、マイナス無限大から指定された...
この関数の差を計算することで、指定範囲の確率（面積）を求...
sageへの入力:
#pre{{
# xが1.2から1.8となる確率は、cum_distribution_functionの...
T.cum_distribution_function(1.8) - T.cum_distribution_fun...
}}
#pre{{
0.0791393511088
}}
*** 正規分布の尤度 [#o03cfe09]
久保本に従って正規分布の尤度をと最尤尤度を求めてみます。
\(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...
& = & \prod_i \frac{1}{\sqrt{2 \pi \sigma^2}} exp\left\{ ...
& = & \prod_i (2 \pi \sigma^2)^{-\frac{1}{2}} exp\left\{ ...
\end{eqnarray}
$$
両辺をlogを取って
$$
\begin{eqnarray}
log L(\mu, \sigma) & = & -\frac{1}{2} \sum_i log(2 \pi \...
& = & -\frac{1}{2} N log(2 \pi \sigma^2) - \frac{1}{2 \si...
\end{eqnarray}
$$
\(\Delta y\)の項を無視すると、
$$
log L(\mu, \sigma) = -\frac{1}{2} N log(2 \pi \sigma^2) -...
$$
となります。これを\(\mu\)で偏微分し、最小となるのは、\(y ...
** ガンマ分布 [#p2cb2fd9]
どんなときにガンマ分布を使うのかGoogleでみると、正の値で...
使われるとありました。
ガンマ分布をSageでプロットしてみます。
sageへの入力:
#pre{{
# Rのdgammaを使用
def _gamma(y, shape, rate):
return N(r('dgamma(%s, %s, %s)' %(y, shape, rate)))
}}
sageへの入力:
#pre{{
g = Graphics()
g += plot(lambda x: _gamma(x, 1, 1), [x, 0, 4], rgbcolor=...
g += plot(lambda x: _gamma(x, 5, 5), [x, 0, 4], rgbcolor=...
#g += plot(lambda x: _gamma(x, 0.5, 0.5), [x, 0, 4], rgbc...
g.show(figsize=5)
}}
&ref(sage4.png);
*** ガンマ分布を使った回帰分析 [#g4e5f04c]
例題では、ある個体の花の重量を\(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 ...
$$
と表され、リンク関数は、expの逆関数のlogとなり、
$$
log \mu_i = a + b log x_i
$$
となります。
ガンマ分布の例題データは、RData形式なので、Rで解析します...
glmを使うと複雑な分布の回帰が簡単に計算できます。しかし、...
プロット結果を見ながら判断するのがよいと実感しました。人...
sageへの入力:
#pre{{
# ガンマ分布用のデータを使って回帰分析
r('load("%s")'%(DATA + 'd.RData'))
r('glm(y ~ log(x), family=Gamma(link="log"), data=d)')
}}
#pre{{
Call: glm(formula = y ~ log(x), family = Gamma(link = "l...
Coefficients:
(Intercept) log(x)
-1.0403 0.6833
Degrees of Freedom: 49 Total (i.e. Null); 48 Residual
Null Deviance: 35.37
Residual Deviance: 17.25 AIC: -110.9
}}
sageへの入力:
#pre{{
# 結果の図化
graph = preGraph("fig-6.13.png")
r('p <- ggplot(data=d, aes(x=x, y=y)) + geom_point() + st...
r('plot(p)')
postGraph(graph)
}}
&ref(fig-6.13.png);
** コメント [#l427012a]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2014/04/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/43/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう：データ解析のための統計モ...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon....
（以下、久保本と書きます）
の第6章の例題をSageを使って再現してみます。
数式処理システムSageのノートブックは、計算結果を表示する...
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#v6dc0e11]
最初に必要なライブラリーやパッケージをロードしておきます。
sageへの入力:
#pre{{
# Rの必要なライブラリ
#r("install.packages('ggplot2')")
r('library(ggplot2)')
r('library(jsonlite)')
# python用のパッケージ
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from ggplot import *
# statsmodelsを使ってglmを計算します
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
** 上限のあるカウントデータの回帰分析 [#e8efce13]
6章の最初に二項分布を使った回帰分析の例がでてきます。ポア...
上限のある場合には、二項分布を使うのだと説明がありました。
（自然界ではこのように上限のあるカウントデータが多いので...
*** 二項分布用のデータについて [#o136f767]
二項分布用のデータは、サポートページにあるdata4a.csvです。
これを読み込んで、データの性質と分布を表示します。
直感的に施肥を施す(f=T)ことでサイズxが大きくなっているの...
また、サイズと種子数の右上がりの関係が見られます。（種子...
sageへの入力:
#pre{{
# 6章のデータを読み込む
d = pd.read_csv(DATA+'data4a.csv')
d.head()
}}
#pre{{
N y x f
0 8 1 9.76 C
1 8 6 10.48 C
2 8 5 10.83 C
3 8 6 10.94 C
4 8 1 9.37 C
[5 rows x 4 columns]
}}
sageへの入力:
#pre{{
d.describe()
}}
#pre{{
N y x
count 100 100.000000 100.000000
mean 8 5.080000 9.967200
std 0 2.743882 1.088954
min 8 0.000000 7.660000
25% 8 3.000000 9.337500
50% 8 6.000000 9.965000
75% 8 8.000000 10.770000
max 8 8.000000 12.440000
[8 rows x 3 columns]
}}
sageへの入力:
#pre{{
ggplot(d, aes(x='x', y='y', color='f')) + geom_point()
}}
#pre{{
<ggplot: (31457393)>
}}
sageへの入力:
#pre{{
ggsave('fig-6.2.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-6.2.png);
*** 二項分布 [#cfca1248]
ここで、Sageのプロット関数を使って二項分布の確率分布がど...
久保本の図6.3と同じ条件でブロットしてみます。
二項分布を
$$
p(y | N, q) = \begin{pmatrix} N \\ y \end{pmatrix} q^y (1...
$$
以下の様に_p関数で定義します。
Sageでプロットしてみます。Sageではグラフ（Graphics）にグ...
とても直感的にグラフを書くことができます。（ggplotも同じ...
sageへの入力:
#pre{{
# 二項分布を定義
def _p(q, y, N):
return binomial(N, y)*q^y*(1-q)^(N-y)
}}
sageへの入力:
#pre{{
g = Graphics()
g += list_plot([_p(0.8, y, 8) for y in (0..8)], plotjoine...
g += list_plot([_p(0.3, y, 8) for y in (0..8)], plotjoin...
g += list_plot([_p(0.1, y, 8) for y in (0..8)], plotjoin...
g.show(figsize=5)
}}
&ref(sage0.png);
*** ロジスティック関数 [#t5a0280e]
ロジスティック関数
$$
q_i = logistic(z_i) = \frac{1}{1 + exp(-z_i)}
$$
の関数を_logisticとして定義し、その分布をSageを使ってプロ...
ロジスティック関数のような曲線を持つデータの場合には、二...
リンク関数は、ロジスティック関数の逆関数であるロジットリ...
ロジット関数は、以下の様に表されます(式a)。
$$
logit(q_i) = log \frac{q_i}{1 - q_i}
$$
sageへの入力:
#pre{{
# ロジスティック関数の定義
def _logistic(z):
return 1/(1 + exp(-z))
}}
sageへの入力:
#pre{{
# ロジスティック曲線
plot(_logistic(x), [x, -6, 6], figsize=5, legend_label='$...
}}
&ref(sage1.png);
*** パラメータ推定 [#c2e896dd]
残念ながら、statmodelsではcbindに相当する二項分布の解析方...
Rを使って計算することにしました。
Sageの中からRの関数を使うことができるので、計算を中断する...
回帰の結果、\(\beta_1 = -19.536, \beta_2 = 1.95, \beta_3=...
sageへの入力:
#pre{{
# statmodelsでcbind(y, N -y)の部分を表現する方法が見つか...
# 0, 1の場合には、statmodelsでも解析可能です。
# 残念ですが、Rを使って処理します。
PandaDf2RDf(d, "d")
}}
sageへの入力:
#pre{{
r('glm(cbind(y, N - y) ~ x + f, data=d, family=binomial)')
}}
#pre{{
Call: glm(formula = cbind(y, N - y) ~ x + f, family = bi...
Coefficients:
(Intercept) x fT
-19.536 1.952 2.022
Degrees of Freedom: 99 Total (i.e. Null); 97 Residual
Null Deviance: 499.2
Residual Deviance: 123 AIC: 272.2
}}
*** 結果の図化 [#z6399d53]
ggplot2のstat_smoothで結果の曲線を表示しようとしたのです...
sageへの入力:
#pre{{
# glmの結果を使って曲線を表示できませんでした！stat_smoot...
graph = preGraph("fig-6.7.png")
#r('p <- ggplot(d, aes(x=x, y=y)) + geom_point() + facet_...
r('p <- ggplot(d, aes(x=x, y=y)) + geom_point() + facet_w...
r('plot(p)')
postGraph(graph)
}}
&ref(fig-6.7.png);
そこで、Sageのプロット機能を使って結果を表示してみました。
sageへの入力:
#pre{{
# 代わりにSageでプロットしてみました。
dT = d[d.f == 'T']
dC = d[d.f == 'C']
dT_pts = list_plot(zip(dT.x, dT.y))
dC_pts = list_plot(zip(dC.x, dC.y), rgbcolor='red')
dT_plt = plot(_logistic(-19.536+1.952*x + 2.022)*8, [x, 7...
dC_plt = plot(_logistic(-19.536+1.952*x)*8, [x, 7, 13], r...
(dT_pts+dT_plt+dC_pts+dC_plt).show(figsize=5)
}}
&ref(sage2.png);
*** 結果の解釈 [#g22c51be]
久保本のすごいところは、単に回帰をすることで終わらず、そ...
リンク関数がロジット関数である(式a)であることから
$$
log \frac{q_i}{1 - q_i} = \beta_1 + \beta_2 x_i + \beta_3...
$$
となり、両辺の指数を取ると、
$$
\frac{q_i}{1 - q_i} = exp(\beta_1 + \beta_2 x_i + \beta_3...
= exp(\beta_1) exp(\beta_2 x_i) exp(\beta_3 f_i)
$$
となります。\(\frac{q_i}{1 - q_i}\)がオッズと呼ばれる量で...
施肥の有無による違いは、exp(2.02)=7.5倍であると推定されま...
** 割り算を使わない方法はすごい！ [#l375b87e]
久保本の6章で感動したのは、割り算を使わない方法です。
私のようなものは、割り算を使って正規化してしまいますが、
このように割り算をしないで、回帰分析する方法としてオフセ...
紹介しています。
オフセット用のデータには、data4b.csvを使います。
このデータは、
- 調査地iの面積を\(A_i\)
- 調査地iの明るさを\(x_i\)
- 調査地iにおける植物個体数\(y_i\)
の記録です。横軸にA, 縦軸にyを取りデータをプロットすると...
面積が大きいほど植物個体数が多く、xが大きい（明るい）程植...
見て取れます。
オフセット項は、offset=np.log(オフセットのカラム)のように...
取った値が使われていますが、これは以下のような理由からで...
人口密度を以下の様に表し、
$$
\frac{平均個体数\lambda_i}{A_i} = 人口密度
$$
平均個体数と明るさに指数関数的な関係があると仮定すると、
$$
\lambda_i = A_i \times 人口密度 = A_i exp(\beta_1 + \beta...
$$
ここで、面積を掛けている部分をlogを使うことで
$$
\lambda_i = exp(\beta_1 + \beta_2 x_i + log A_i)
$$
と変形することができます。つまり、\(log A_i\)だけずれた値...
よいとありました。
sageへの入力:
#pre{{
# オフセット用のデータを読み込む
d2 = pd.read_csv(DATA+'data4b.csv')
d2.head()
}}
#pre{{
y x A
0 57 0.68 10.3
1 64 0.27 15.6
2 49 0.46 10.0
3 64 0.45 14.9
4 82 0.74 14.0
[5 rows x 3 columns]
}}
sageへの入力:
#pre{{
ggplot(d2, aes(x='A', y='y', color='x')) + geom_point()
}}
#pre{{
<ggplot: (40327389)>
}}
sageへの入力:
#pre{{
ggsave('fig-6.10-A.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-6.10-A.png);
*** 推定値 [#b4892baf]
推定された値をプロットすると以下のようになります。
データよりもはっきりと明るさによる影響がきれいに分かりま...
sageへの入力:
#pre{{
fit = smf.glm('y ~ x', data=d2, offset=np.log(d.A), famil...
fit.summary2()
d2['pred'] = fit.predict()
}}
sageへの入力:
#pre{{
# 予測値のプロット
ggplot(d2, aes(x='A', y='pred', color='x')) + geom_point()
}}
#pre{{
<ggplot: (39619317)>
}}
sageへの入力:
#pre{{
ggsave('fig-6.10-C.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-6.10-C.png);
*** オフセット項を使わないと [#k92b37cd]
オフセット項を使わないと面積による影響が回帰分析に含まれ...
sageへの入力:
#pre{{
# 人口密度Aをオフセットを使わない場合
fit_2 = smf.glm('y ~ x', data=d2, family=sm.families.Pois...
fit_2.summary2()
d2['pred2'] = fit_2.predict()
}}
sageへの入力:
#pre{{
ggplot(d2, aes(x='A', y='pred2', color='x')) + geom_point()
}}
#pre{{
<ggplot: (40798181)>
}}
sageへの入力:
#pre{{
ggsave('fig-6.10-D.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-6.10-D.png);
** 正規分布と尤度 [#m8c6dfd3]
Sageにも確率分布を扱う関数RealDistirubtionがVer.5から導入...
これを使って、正規分布の確率密度関数をプロットしていまし...
sageへの入力:
#pre{{
# Sage ver.5から導入されたRealDistributionを使う
T1 = RealDistribution('gaussian', 1, seed=100)
T2 = RealDistribution('gaussian', 3, seed=100)
}}
sageへの入力:
#pre{{
T1_plt = plot(lambda x : T1.distribution_function(x), [x,...
T2_plt = plot(lambda x : T2.distribution_function(x), [x,...
T3_plt = plot(lambda x : T1.distribution_function(x-2), [...
(T1_plt + T2_plt + T3_plt).show(figsize=5)
}}
&ref(sage3.png);
*** 確率密度関数から確率を得る [#u6ff9651]
確率密度関数から指定された範囲に含まれる確率は、その範囲...
cum_distribution_functionは、マイナス無限大から指定された...
この関数の差を計算することで、指定範囲の確率（面積）を求...
sageへの入力:
#pre{{
# xが1.2から1.8となる確率は、cum_distribution_functionの...
T.cum_distribution_function(1.8) - T.cum_distribution_fun...
}}
#pre{{
0.0791393511088
}}
*** 正規分布の尤度 [#o03cfe09]
久保本に従って正規分布の尤度をと最尤尤度を求めてみます。
\(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...
& = & \prod_i \frac{1}{\sqrt{2 \pi \sigma^2}} exp\left\{ ...
& = & \prod_i (2 \pi \sigma^2)^{-\frac{1}{2}} exp\left\{ ...
\end{eqnarray}
$$
両辺をlogを取って
$$
\begin{eqnarray}
log L(\mu, \sigma) & = & -\frac{1}{2} \sum_i log(2 \pi \...
& = & -\frac{1}{2} N log(2 \pi \sigma^2) - \frac{1}{2 \si...
\end{eqnarray}
$$
\(\Delta y\)の項を無視すると、
$$
log L(\mu, \sigma) = -\frac{1}{2} N log(2 \pi \sigma^2) -...
$$
となります。これを\(\mu\)で偏微分し、最小となるのは、\(y ...
** ガンマ分布 [#p2cb2fd9]
どんなときにガンマ分布を使うのかGoogleでみると、正の値で...
使われるとありました。
ガンマ分布をSageでプロットしてみます。
sageへの入力:
#pre{{
# Rのdgammaを使用
def _gamma(y, shape, rate):
return N(r('dgamma(%s, %s, %s)' %(y, shape, rate)))
}}
sageへの入力:
#pre{{
g = Graphics()
g += plot(lambda x: _gamma(x, 1, 1), [x, 0, 4], rgbcolor=...
g += plot(lambda x: _gamma(x, 5, 5), [x, 0, 4], rgbcolor=...
#g += plot(lambda x: _gamma(x, 0.5, 0.5), [x, 0, 4], rgbc...
g.show(figsize=5)
}}
&ref(sage4.png);
*** ガンマ分布を使った回帰分析 [#g4e5f04c]
例題では、ある個体の花の重量を\(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 ...
$$
と表され、リンク関数は、expの逆関数のlogとなり、
$$
log \mu_i = a + b log x_i
$$
となります。
ガンマ分布の例題データは、RData形式なので、Rで解析します...
glmを使うと複雑な分布の回帰が簡単に計算できます。しかし、...
プロット結果を見ながら判断するのがよいと実感しました。人...
sageへの入力:
#pre{{
# ガンマ分布用のデータを使って回帰分析
r('load("%s")'%(DATA + 'd.RData'))
r('glm(y ~ log(x), family=Gamma(link="log"), data=d)')
}}
#pre{{
Call: glm(formula = y ~ log(x), family = Gamma(link = "l...
Coefficients:
(Intercept) log(x)
-1.0403 0.6833
Degrees of Freedom: 49 Total (i.e. Null); 48 Residual
Null Deviance: 35.37
Residual Deviance: 17.25 AIC: -110.9
}}
sageへの入力:
#pre{{
# 結果の図化
graph = preGraph("fig-6.13.png")
r('p <- ggplot(data=d, aes(x=x, y=y)) + geom_point() + st...
r('plot(p)')
postGraph(graph)
}}
&ref(fig-6.13.png);
** コメント [#l427012a]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
