Dobson(序論)
Yuta NAGANO
11/1/2019
※ 全体的に重要なことは書いていないので、軽く目を通す程度で良いと思います。
1.1 背景
- 特に重要なことは書いてなかった。
1.2 範囲
- Dobson本で取り上げられる統計手法は、様々な測定値間の関係を分析するための手法
説明変数（予測・独立変数）
- 実験計画によって固定された値のようにランダムでない測定値または観測値としてあつかう。
応答変数（結果・従属変数）
- 説目変数の値に応じて変動する測定値であり、確率変数としてあつかう。
変数の尺度
説明変数および応答変数は、以下のような尺度のいずれかを用いて測定される。
名義尺度
- 順序のないカテゴリー
例）オス・メスや生存・死亡など
順序尺度
- 何らかの自然的な順序が存在するカテゴリー
例）小学生の学年
連続的な測定値
- 少なくとも理論的に、ある連続体上のどこにでも落ちる可能性のある値
例）重さや時間、比率
尺度の上位グループ
名義・順序尺度はカテゴリー変数や離散変数と呼び、カテゴリーに対する観測値の数、つまりカウントデータや頻度データが記録される。
連続的な測定値については連続変数と呼び、ここの連続データが記録される。
特に、連続変数のことを量的、離散変数のことを質的と呼ぶことがしばしばある。
質的な説明変数を因子（factor）、そのカテゴリーを水準（level）と呼び、量的な説明変数を共変量（covariate）と呼ぶ。
1.3 表記法
一般的に、確率変数をイタリックのアルファベット大文字、観測値を対応する小文字で表す。
- 確率変数：\(Y_1, Y_2, Y_3, ..., Y_n\)
- 観測値：\(y_1, y_2, y_3, .... y_n\)
また、ギリシャ文字はパラメータを、それに対応するイタリックのアルファベットの小文字は推定量または推定値を表す。
場合により（こっちの方が多い気がしますが…）、記号\(\hat{}\)を推定量または推定値の表現に使われる。
- 推定したいパラメータ：\(\beta{}\)
- パラメータの推定値（推定量）：\(b\)または\(\hat{\beta}\)
ベクトルと行列は、ランダムであるかどうかに関わらず、それぞれ太文字を使う。
- 観測or確率変数ベクトル：\(\boldsymbol{y} = \begin{bmatrix}y_1\\.\\.\\.\\y_n\end{bmatrix} or \begin{bmatrix}Y_1\\.\\.\\.\\Y_n\end{bmatrix}\)
- パラメータベクトル：\(\boldsymbol{\beta} = \begin{bmatrix}\beta_1\\.\\.\\.\\\beta_n\end{bmatrix}\)
- 行列：\(\boldsymbol{X}\)
また、上付き添字の\(^T\)は行列の転置を表し、列ベクトルを行ベクトルで表したい時は、\(\boldsymbol{y} = [Y_1, ..., Y_n]^T\)とする。
連続的な確率変数\(Y\)の確率密度関数（\(Y\)が離散変数の場合、確率関数）を以下のように表記する。
\[f(y; \theta{})\] ここで、\(\theta\)は分布パラメータであり、必ずしも１つとは限らない。 また、和を表すために\(.\)、平均を表すために\(\bar{}\)を用いて、以下のように表現できる。
\[\bar{y} = \frac{1}{N}\sum^N_{i = 1}{y_i} = \frac{1}{N}y.\]
確率変数\(Y\)の期待値（平均値）と分散はそれぞれ、\(E(Y)\)、\(var(Y)\)と表される。
1.4 正規分布および正規分布から導かれる分布
多くの推定量や検定統計量の標本分布は正規分布と関係がある。
それらは、正規分布に従う確率変数として関係するか、中心極限定理を通して漸近的に関係するか。
1.4.1 正規分布
確率変数\(Y\)が平均\(\mu{}\)、分散\(\sigma^2\)の正規分布に従うとき、その確率密度関数は、 \[f(y; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}exp[-\frac{1}{2}(\frac{y-\mu}{\sigma})^2]\] であり、\(Y \sim N(\mu, \sigma^2)\)と表す。
- 平均0、分散1の正規分布を標準正規分布と呼ぶ（\(Y \sim N(0, 1)\)）。
- \(Y_1, ..., Y_n\)を\(Y_i \sim N(\mu_i, \sigma_i^2), i = 1, ..., n\)であるような確率変数、\(Y_i\)と\(Y_j\)の共分散を
\[cov(Y_i, Y_j) = \rho_{ij}\sigma_i\sigma_j\]
と表す。ここで、\(\rho_{ij}\)は\(Y_i\)と\(Y_j\)の相関係数を表す。
この時、\(Y_1, ..., Y_n\)の同時分布は平均ベクトルとして、\(\boldsymbol{\mu} = [\mu_1, ..., \mu_n]^T\)、共分散行列として対角要素が\(\sigma^2_i\)、非対角要素が\(\rho_{ij}\sigma_i\sigma_j(i \neq j)\)であるような行列\(\boldsymbol{V}\)をもつ多変量正規分布である。
このことを、\(\boldsymbol{y} \sim \boldsymbol{N}(\boldsymbol{\mu}, \boldsymbol{V})\)とかき、\(\boldsymbol{y} = [Y_1, ..., Y_n]^T\)である。
確率変数\(Y_1, ..., Y_n\)は独立に正規分布に従い、\(Y_i \sim N(\mu_i, \sigma_i^2), i = 1, ..., n\)であり、 \[W = a_1Y_1+a_2Y_2+...+a_nY_n\] とする（\(a_1, ..., a_n\)は定数）。この時、\(W\)も正規分布に従い（正規分布の和は正規分布）、 \[W = \sum^n_{i=1}{a_iY_i} \sim N(\sum^n_{i=1}a_i\mu_i, \sum^b_{i=1}a_i^2\sigma_i^2)\] である。
２つの正規乱数（正規分布）の和は正規分布なのか with R
set.seed(123)
# 正規乱数を発生させるために平均と標準偏差を定義
mu1 <- 2
sd1 <- 1
mu2 <- -1
sd2 <- 0.4
# それぞれの正規分布からのランダムサンプリング
(data <- tibble(d_norm1 = rnorm(1000, mu1, sd1),
d_norm2 = rnorm(1000, mu2, sd2)))## # A tibble: 1,000 x 2
## d_norm1 d_norm2
## <dbl> <dbl>
## 1 1.44 -1.40
## 2 1.77 -1.42
## 3 3.56 -1.01
## 4 2.07 -1.05
## 5 2.13 -2.02
## 6 3.72 -0.584
## 7 2.46 -0.900
## 8 0.735 -0.0335
## 9 1.31 -0.726
## 10 1.55 -1.18
## # … with 990 more rows# それぞれの正規乱数をヒストグラムとして描画
ggplot(data) + geom_histogram(aes(x = d_norm1), fill = "#F8766D", color = "#F8766D", alpha = 0.5) +
geom_histogram(aes(x = d_norm2), fill = "#00BFC4", color = "#00BFC4", alpha = 0.5)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# それぞれの和のヒストグラムを描画
ggplot(data) +
geom_histogram(aes(x = d_norm1 + d_norm2), fill = "#7CAE00", color = "#7CAE00", alpha = 0.5)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# 最後に和の正規性の検定
shapiro.test(x = (data$d_norm1 + data$d_norm2))##
## Shapiro-Wilk normality test
##
## data: (data$d_norm1 + data$d_norm2)
## W = 0.99859, p-value = 0.61641.4.2 カイ２乗分布
- 自由度nの中心カイ二乗分布は、標準正規分に従うn個の独立な確率変数\(Z_i, ..., Z_n\)の２乗和の分布として定義される。このことを、 \[X^2 = \sum^n_{i=1}Z^2_i \sim \chi^2(n)\] のように表す。
行列を用いた表記法では、\(\boldsymbol{z}=[Z_i, ..., Z_n]^T\)ならば\(\boldsymbol{z}^T\boldsymbol{z} = \sum^n_{i=1}Z_i^2\)なので、\(X^2 = \boldsymbol{z}^T\boldsymbol{z}\sim\chi^2(n)\)と表す。
- \(X^2\)が自由度nのカイ２乗分布\(\chi^2(n)\)に従う時、期待値は\(E(X^2) = n\)、分散は\(var(X^2) = 2n\)である。
- \(Y_1, ..., Y_n\)が独立に正規分布に従う確率変数で、\(Y_i \sim N(\mu_i, \sigma_i^2)\)である時、 \[X^2 = \sum^n_{i = 1}(\frac{Y_i - \mu_i}{\sigma_i})^2 \sim \chi^2(n)\] である。（変数\(Z_i = (Y_i-\mu_i)/\sigma_i\)が\(N(0,1)\)に従うから）
- \(Z_1, ..., Z_n\)が分布\(N(0,1)\)に従う独立に確率変数で、\(Y_i = Z_i + \mu_i\)とする。
ただし、\(\mu_i, ..., \mu_n\)は少なくとも１つがゼロ出ない。
この時、 \[\sum Y_i^2 = \sum (Z_i + \mu_i)^2 = \sum Z_i^2 + 2\sum Z_i\mu_i + \sum \mu_i^2\] の分布は\(\chi^2(n)\)に比べて、平均も分散も大きく、それぞれ\(n + \lambda\)、\(2n + 4\lambda\)である。
ただし、\(\lambda = \sum \mu_i^2\)とする。この分布を自由度n、非心パラメータ\(\lambda\)の非心カイ２乗分布と呼び、\(\chi^2(n, \lambda)\)と表わす。
- \(Y_i, ..., Y_n\)が必ずしも互いに独立ではなく、ベクトル\(\boldsymbol{y} \sim \boldsymbol{N}(\boldsymbol{\mu},\boldsymbol{V})\)とする。
ただし、共分散行列\(\boldsymbol{V}\)は正則で、逆行列\(\boldsymbol{V}^{-1}\)をもつ。
この時、 \[X^2 = (\boldsymbol{y}-\boldsymbol{\mu})^T\boldsymbol{V}^{-1}(\boldsymbol{y}-\boldsymbol{\mu}) \sim \chi^2(n)\]である。 - 一般に、\(\boldsymbol{y} \sim \boldsymbol{N}(\boldsymbol{\mu},\boldsymbol{V})\)である時、確率変数\(\boldsymbol{y}^T\boldsymbol{V}^{-1}\boldsymbol{y}\)は非心カイ２乗分布\(\chi^2(n, \lambda)\)に従う。ただし、\(\lambda = \boldsymbol{\mu}^T\boldsymbol{V}^{-1}\boldsymbol{\mu}\)である。
- \(X^2_1, ..., X^2_n\)が中心または非心カイ２乗分布に従う互いに独立な確率変数で、\(X_i^2 \sim \chi^2(n_i, \lambda_i)\)である時、それらの和は自由度\(\sum n_i\)、非心パラメータ\(\sum \lambda_i\)のカイ２乗分布に従う。 すなわち、 \[\sum^m_{i = 1} X^2_i \sim \chi^2(\sum^m_{i = 1} n_i, \sum^m_{i=1}\lambda_i)\] である。この性質は、カイ２乗分布の再生性と呼ばれる。
- \(\boldsymbol{y} \sim \boldsymbol{N}(\boldsymbol{\mu},\boldsymbol{V})\), \(\boldsymbol{y}\)はn個の要素をもち、\(\boldsymbol{V}\)は階数\(k < n\)で特異、
したがって、逆行列が一意には定まらないものとする。
このとき、\(\boldsymbol{V}\)の一般逆行列を\(\boldsymbol{V}^-\)と表わすと、確率変数\(\boldsymbol{y}^T\boldsymbol{V}^{-}\boldsymbol{y}\)は自由度\(k\)、非心パラメータ\(\lambda = \boldsymbol{\mu}^T\boldsymbol{V}^{-1}\boldsymbol{\mu}\)の非心カイ２乗分布に従う。
1.4.3 t分布
- 自由度nのt分布は、２つの独立な確率変数の比として定義される。
分子は標準正規分布に従う確率変数、分母は中心カイ２乗分布に従う確率変数をその自由度で割って平方根をとったものである。
つまり、 \[T = \frac{Z}{(X^2/n)^{1/2}}\] である。
ここに\(Z \sim N(0, 1), X^2 \sim \chi^2(n)\)であり、\(Z\)と\(X^2\)は独立とする。
このことを$ T t(n)$
1.4.4 F分布
- 自由度nとmをもつ中心F分布は２つの独立な中心カイ２乗確率変数をそれぞれの自由度で割ったものの比として定義される。 \[F = {\frac{X^2_1}{n}}/{\frac{X^2_2}{m}}\] ただし、\(X_1^2 \sim \chi^2(n), X_2^2 \sim \chi^2(m)\)であり、\(X_1^2\)と\(X^2_2\)は互いに独立であるとする。
上のことを\(F ~ F(n, m)\)と表わす。
- t分布とF分布の関係として、t分布の式を２乗し、上記式を用いることで、 \[T^2 = \frac{Z^2}{1}/\frac{X^2}{n}\sim F(1,n)\] が得られる。
すなわち自由度nのt分布\(t(n)\)に従う確率変数をそれぞれの自由度で割ったものの比として定義される。
分子は非心カイ２乗分布、分母は中心カイ２乗分布である。
すなわち \[F = \frac{X^2_1}{n}/\frac{X^2_2}{m}\] ここで、\(X^2_1 \sim \chi^2(n, \lambda), X^2_2 \sim \chi^2(m)\)であり、\(X^2_1\)と\(X^2_2\)は互いに独立とする。
非心F分布の平均は同じ自由度をもつ中心F分布の平均よりも大きくなる。
1.5 二次形式
- 二次形式とは、全ての項の次数が２である多項式。
例）\(y^2_1+y_2^2 + y_3^2\) - \(\boldsymbol{A}\)を次のような対象行列とおく、すなわち、\(a_{ij} = a{ji}\)とする。 \[\begin{bmatrix}a_{11}&a_{12}&.&.&.&a_{1n}\\a_{21}&a_{22}&.&.&.&a_{2n}\\.&&.&&&.\\.&&&.&&.\\.&&&&.&.\\a_{n1}&a_{n2}&.&.&.&a_{nn} \end{bmatrix}\] この時、式\(\boldsymbol{y}^T\boldsymbol{Ay} = \sum_i\sum_j a_{ij}y_iy_j\)は\(y_i\)の二次形式である。
\((\boldsymbol{y} - \boldsymbol{\mu})^T\boldsymbol{V}^{-1}(\boldsymbol{y}-\boldsymbol{\mu})\)は\((y_i-\mu_i)\)の二次形式であるが、\(y_i\)の二次形式ではない。
- 二次形式\(\boldsymbol{y}^T\boldsymbol{Ay} = \sum_i\sum_j a_{ij}y_iy_j\)と行列\(\boldsymbol{A}\)は、要素の全てが0でない任意の\(\boldsymbol{y}\)に対して、\(\boldsymbol{y}^T\boldsymbol{Ay} > 0\)である時、正定値であると言われる。
正定値であるための必要十分条件は次の全ての行列式が正となることである。 \[|A_1| = a_{11}, |A_2| = \begin{vmatrix} a_{11}&a_{12}\\a_{21}&a_{22}\end{vmatrix}, |A_3| = \begin{vmatrix} a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{vmatrix}, |A_n| = det\boldsymbol{A}\]
- 行列\(\boldsymbol{A}\)の階数は二次形式\(Q = \boldsymbol{y}^T\boldsymbol{Ay}\)の自由度とも呼ばれる。
- \(Y_1, ..., Y_n\)をそれぞれ正規分布\(N(0, \sigma^2)\)に従う独立な確率変数と仮定する。
\(Q = \sum^n_{i=1} Y_i^2\)と起き、\(Q_1, ..., Q_k\)を\(Y_i\)の二次形式として次の関係が成り立つとする。
\[Q = Q_1 + Q_2 + ... + Q_k\] ここで\(Q_i\)は自由度\(m_i (i = 1, ..., k)\)をもつ。
すると\(m_1 + m_2 + ... + m_k = n\)の時、かつその時に限り、\(Q_1, ..., Q_k\)は独立な確率変数であり、\(Q_1/\sigma^2 \sim \chi^2(m_1), Q_2/\sigma^2 \sim \chi^2(m_2), ..., Q_k/\sigma_k \sim \chi^2(m_k)\)となる。
これはCochranの定理である。
- Cochranの定理より、次が成り立つ。２つの確率変数\(X^2_1 \sim \chi^2(m), X^2_2 \sim\chi^2(k)\)(ただし、\(m > k\))があるとする。その差が非負定値、すなわち\(X^2 = X_1^2-X_2^2 \geq 0\)ならば、\(X^2\sim \chi^2(m - k)\)である。
1.6 推定
1.6.1 最尤推定
- \(\boldsymbol{y} = [Y_1, ..., Y_n]^T\)を確率変数ベクトルとし、\(Y_i\)の同時確率密度関数を次のように表わす。
\[f(y; \boldsymbol{\theta})\] ただし、\(\boldsymbol{\theta}\)はパラメータベクトル\(\boldsymbol{\theta} = [\theta_1, ..., \theta_p]^T\)である。
- 尤度関数\(L(\boldsymbol{\theta}; y)\)は代数的には同時確率密度関数\(f(y; \boldsymbol{\theta})\)と同じであるが、表記が\(\boldsymbol{\theta}\)を固定した時の確率変数\(\boldsymbol{y}\)を考えることから、\(\boldsymbol{y}\)を固定した時のパラメータ\(\boldsymbol{\theta}\)を考えるように変更されている。\(L\)は確率変数ベクトル\(\boldsymbol{y}\)の関数として定義され、それ自体が確率変数である。
- \(\Omega\)をパラメータベクトル\(\boldsymbol{\theta}\)の可能な全ての値を含む集合と定義し、\(\Omega\)をパラメータ空間と呼ぶ。
- \(\boldsymbol{\theta}\)の最尤推定量は尤度関数を最大にする\(\boldsymbol{\hat{\theta}}\)である。
すなわち、 \[L(\boldsymbol{\hat{\theta}}; y) \geq L(\boldsymbol{\theta}; y),\;\;\boldsymbol{\theta}\in\Omega\] - 対数関数は単調関数であるので、\(\boldsymbol{\hat{\theta}}\)は対数尤度関数\(l(\boldsymbol{\theta}; \boldsymbol{y}) = logL(\boldsymbol{\theta}; \boldsymbol{y})\)を最大にする値でもある。 \[l(\boldsymbol{\hat{\theta}}; \boldsymbol{y}) \geq l(\boldsymbol{\theta}; \boldsymbol{y}), \;\; \boldsymbol{\theta} \in \Omega\] 基本的には、尤度関数よりも対数尤度関数の方が計算が楽。
- 通常、推定量\(\boldsymbol{\hat{\theta}}\)は対数尤度関数を\(\boldsymbol{\theta}\)の各要素\(\boldsymbol{\theta}_j\)に関して、偏微分して0とおいた、以下の連立方程式をとくことにより得られる。 \[\frac{\partial{l}(\boldsymbol{\theta}; \boldsymbol{y})}{\partial\theta_j} = 0, \;\;j = 1, ..., p\] そして、上記式の解が\(l(\boldsymbol{\theta}; \boldsymbol{y})\)の最大値に対応するかの確認には、２階導関数 \[\frac{\partial^2l(\boldsymbol{\theta}; \boldsymbol{y})}{\partial\boldsymbol{\theta}_j\partial\boldsymbol{\theta}_k}\] を\(\boldsymbol{\theta} = \boldsymbol{\hat{\theta}}\)とおいて計算した行列が負定値行列となっているかを調べる必要がある。
\(\boldsymbol{\theta}\)が１次元で１つだけの要素\(\theta\)をもつ場合、以下の式を確認する必要がある。 \[[\frac{\partial^2l(\boldsymbol{\theta}; \boldsymbol{y})}{\partial\theta^2}]_{\theta = \hat{\theta}} < 0\] - パラメータ空間\(\Omega\)の叶華状の値\(\theta\)で\(l(\boldsymbol{\theta}; \boldsymbol{y})\)の局所最大値（極値）を与える\(\boldsymbol{\theta}\)があるかどうかの確認も必要である。
全ての局所最大値を調べてその中で最大となる値\(\boldsymbol{\hat{\theta}}\)が最尤推定量となる。 - 最尤推定量の重要な性質の１つに\(g(\boldsymbol{\theta})\)をパラメータ\(\boldsymbol{\theta}\)の任意の関数とする時、\(g(\boldsymbol{\theta})\)の最尤推定量は\(g(\boldsymbol{\hat{\theta}})\)となるという性質があり、最尤推定量の不変性と呼ばれる。
つまり、最尤推定を行うのに便利なパラメータの関数に対して尤度関数を最大化し、最尤推定量の不変性を利用して必要なパラメータの最尤推定量を求めれば良い。
1.6.2 例：ポアソン分布
- \(Y_1, ..., Y_n\)を以下に示すような同じパラメータ\(\theta\)をもつポアソン分布に従う独立な確率変数とする。 \[f(y_i; \theta) = \frac{\theta^{y_i}e^{-\theta}}{y!}, \;\; y_i = 0, 1, 2, ...\]
- その同時分布は \[\begin{align} f(y_1, ..., y_n; \theta) &= \prod^n_{i = 1}f(y_i; \theta) = \frac{\theta^{y_1}e^{-\theta}}{y_1!}×\frac{\theta^{y_2}e^{-\theta}}{y_2!}×...×\frac{\theta^{y_n}e^{-\theta}}{y_n!} \\ &= \frac{\theta^{\sum y_i}e^{-n\theta}}{y_1!y_2!...y_n!} \end{align}\] のようになり、尤度関数\(L(\theta; y_1, ..., y_n)\)もまた同じ式で表される。
最尤推定量\(\hat{\theta}\)を見つけるためには対数尤度関数 \[l(\theta; y_1, ..., y_n) = logL(\theta; y_1, ..., y_n) = (\sum y_i)log\theta-n\theta-\sum(logy_i!)\] を利用する方が便利である。\(l\)の導関数は、 \[\frac{dl}{d\theta} = \frac{1}{\theta}\sum{\frac{y_i}{n}} = \bar{y}\]
ポアソン乱数のパラメータ最尤推定 with R
set.seed(123)
# ポアソン乱数を発生させる
L <- 3.4
d_pois <- rpois(1000, L)
# 負の対数尤度を計算する
logL <- function(lambda) -sum(dpois(d_pois, lambda, log = TRUE))
# 負の対数尤度を最適化する
stats4::mle(logL, start = list(lambda = 1))##
## Call:
## stats4::mle(minuslogl = logL, start = list(lambda = 1))
##
## Coefficients:
## lambda
## 3.392# データの平均を計算
mean(d_pois)## [1] 3.3921.6.3 最小２乗法
- \(Y_1, ..., Y_n\)をそれぞれ期待値\(\mu_1, ..., \mu_n\)をもつ独立な確率変数とし、各期待値はパラメータベクトル\(\boldsymbol{\beta} = [\beta_1, ..., \beta_n]^T, p < n\)の関数でありとする。
すなわち、
\[E(Y_i) = \mu_i(\boldsymbol{\beta})\] 最小二乗法のもっとも簡単な形は、観測値\(Y_i\)と期待値\(\mu_i\)の差の平方和を最小にする推定量\(\hat{\beta}\)を見つけることである。
\[S = \sum [Y_i - \mu_i(\beta)]^2\] 通常\(\hat{\beta}\)は\(S\)はを各要素\(\beta_j\)に関して微分し、連立方程式
\[\frac{\partial S}{\partial \beta_j} = 0, \;\; j = 1, ..., p\] を解くことにより得られる。
- もちろん解が最小値に対応しているかを確認する必要がある（２回微分した行列が正定値行列であるか）。これらの解とパラメータ空間の境界での局所最小値の中から大域的最小値を求める必要がある。
- ここで\(Y_i\)はそれぞれ全てが等しくはない分散\(\sigma^2\)をもつとしよう。このとき以下の重み付き平方和を最小にする方が望ましい。
\[S = \sum w_i[Y_i - \mu_i(\boldsymbol{\beta})]^2\] ただし、重みは\(w_i = (\sigma_i^2)^{-1}\)である。このようにすると、信頼性が低い観測値（\(Y_i\)が大きい分散をもつ観測値）は推定値に対して影響が少なくなる。
- より一般的に、\(\boldsymbol{y} = [Y_1, ..., Y_n]^T\)を平均ベクトル\(\boldsymbol{\mu} = [\mu_1, ..., \mu_n]^T\)と分散共分散行列\(\boldsymbol{V}\)をもつ確率変数ベクトルとする。重み付き最小２乗推定量は以下のSを最小化することによって得られる。
\[S = (\boldsymbol{y}-\boldsymbol{\mu})^T\boldsymbol{V}^{-1}(\boldsymbol{y}-\boldsymbol{\mu})\] - 本書本章に書いていないが、最小二乗法は正規分布に従う変数にしか適用できない。また、正規分布に従う変数のパラメータを最尤推定した場合、最小二乗法と同じ値になることがわかっている（理論的にも）。
1.6.4 推定に関する注釈
特になし
1.6.5 例：サイクロン（熱帯低気圧）with R
データの作成
(data <- tibble(season = seq(1:13),
cyclone = c(6, 5, 4, 6, 6, 3, 12, 7, 4, 2, 6, 7, 4)))## # A tibble: 13 x 2
## season cyclone
## <int> <dbl>
## 1 1 6
## 2 2 5
## 3 3 4
## 4 4 6
## 5 5 6
## 6 6 3
## 7 7 12
## 8 8 7
## 9 9 4
## 10 10 2
## 11 11 6
## 12 12 7
## 13 13 4サイクロン数は0以上無限大の離散値であるため、ポアソン分布に従う確率変数であると言える。
また、ポアソン分布の唯一のパラメータlambdaは、平均値である。
mean(data$cyclone)## [1] 5.538462もう１つのアプローチとして、最尤推定がある。
logL <- function(lambda) -sum(dpois(data$cyclone, lambda, log = TRUE))
stats4::mle(logL, start = list(lambda = 1))## Warning in dpois(data$cyclone, lambda, log = TRUE): 計算結果が NaN になりま
## した##
## Call:
## stats4::mle(minuslogl = logL, start = list(lambda = 1))
##
## Coefficients:
## lambda
## 5.538548lambdaが最尤推定値のポアソン分布を描画する。
set.seed(123)
tibble(x = c(1:25), y = dpois(x, mean(data$cyclone))) %>% ggplot() +
geom_line(aes(x = x, y = y)) + theme_bw()