Def: 代表値(Averages)
▶ 標本平均
データ \(X = \{x_1, \cdots, x_n\}\) が与えられとき，標本平均（sample mean） \(\overline{x}\) は次になります：
\[ \overline{x} = \frac{x_1 + \cdots + x_n}{n} \]
標本平均は分布の代表値として最も使用されるものだが，外れ値に対して弱い性質がある．
import numpy as np
X_0 = np.array([5.6, 5.7, 5.4, 5.5, 5.8, 5.2, 5.3, 5.6, 5.4, 55.5])
X_1 = np.array([5.6, 5.7, 5.4, 5.5, 5.8, 5.2, 5.3, 5.6, 5.4])
print(
"""X_0: sample mean = {}, median = {}\nX_1: sample mean = {}, median = {}
""".format(
np.mean(X_0), np.median(X_0), np.mean(X_1), np.median(X_1)
)
)X_0: sample mean = 10.5, median = 5.55
X_1: sample mean = 5.5, median = 5.5
上記の例のように，medianは外れ値の混入があってもその影響は軽微ですが，標本平均は大きく変わっており外れ値に対して弱いことがわかる．
▶ 刈り込み平均
外れ値の影響を弱めて標本平均を推定する方法として，刈り込み平均(trimmed mean)があります． 上側 \(100\alpha \%\) と下側 \(100\alpha \%\) を使わないで推定する方法で，\(x_{[i]}\) を順序統計値として
\[ \hat\mu_\alpha = \frac{1}{n-2m} \sum_{i=m+1}^{n-m}x_{[i]}, \ \ m = \lfloor n\alpha \rfloor \]
で推定する方法を刈り込み平均という．利用にあたって，外れ地の割合を事前に想定する必要がありますが， 少々適当に推定しても妥当な推定になりやすい特徴があります．
\(\alpha = 0.1\) としてPythonで計算してみると以下，
from scipy import stats
self_trimmed_mean = np.mean([5.7, 5.4, 5.5, 5.8, 5.4, 5.3, 5.6, 5.6])
trimmed_mean = stats.trim_mean(X_0, 0.1)
print(self_trimmed_mean, trimmed_mean)5.5375 5.5375Trimmed meanはARE(Asymptotic relative efficiency, 漸近相対効率)という観点からも大抵の裾の重さに対して（重すぎるのは厳しいですが．．．）高いパフォーマンスがあることが知られています．
▶ 幾何平均
\(x_i > 0\) となるようなデータについて，幾何平均は以下のように計算されます：
\[ \overline{x}_G = \bigg(\prod_{i=1}^n x_i \bigg)^{\frac{1}{n}} \]
2000年から2005年までのXこくのでの物価上昇率が2%, 5%, 2%, 5%, 10%とあるとき，年平均上昇率は算術平均ではなく幾何平均で計算すべきで
\[ (1.02 \times 1.05 \times 1.02 \times 1.05 \times 1.1)^{1/5} \approx 1.0476 \]
すなわち年平均約4.8%の増加と報告すべきとなります．なお，相加相乗平均より，幾何平均は算術平均より小さい値になることがわかります．もし大きい値を出してしまっていたら計算ミスを疑うべきです．
幾何平均について対数をとると以下のように算術平均で表すことができます
\[ \log(\overline{x}_G) = \frac{1}{n}\sum \log(x_i) \]
ここから，幾何平均を計算するときは一旦log transformationを実行し，算術平均を計算し，その後 \(\exp(\cdot)\) で元のスケールに戻すという形でよく計算されます．
クラス分類の評価指標との関係では，sensitivity(感応度)とspecificity(特異度)の幾何平均を用いたG-Mean(geometric mean)という指標があります．
\[ \begin{align*} \operatorname{G-mean} &= \sqrt{\operatorname{sensitivity} \times \operatorname{specificity}}\\ &= \sqrt{\operatorname{recall} \times \operatorname{True Negative Rate}} \end{align*} \]
▶ 調和平均
\(x_i > 0\) となるようなデータについて，調和平均(harmonic mean)は以下のように計算されます：
\[ \frac{1}{\overline{x}_H} = \frac{1}{n}\sum\frac{1}{x_i} \]
とある車が距離 \(\alpha\) の区間Aでは25km/h, 距離 \(\beta\) の区間Bでは15km/hで走っていたとします．このとき，この車の平均時速は
\[ \frac{1}{\text{平均時速}} = \frac{\alpha}{\alpha + \beta} \frac{1}{25} + \frac{\beta}{\alpha + \beta} \frac{1}{15} \]
\(\alpha = \beta\) のときは
\[ \frac{1}{\text{平均時速}} = \frac{1}{2} \bigg(\frac{1}{25} + \frac{1}{15}\bigg) \]
平均を計算するにあたって，値が同じスケールの単位である必要であるため，上の平均時速の例では調和平均を利用することが 好ましいとされます．なお区間Aをx時間で25km/h, 区間Bをy時間で15km/hという場合はウェイトが時間単位で表されているので
\[ \text{平均時速} = \frac{x}{x + y} \times 25 + \frac{y}{x + y} \times 15 \]
モデルの評価指標の１つにprecisionとrecallを用いたF1-scoreがありますが，precisionとrecallも分子はそれぞれTrue Positiveで共通していますが，分母がそれぞれ \(\operatorname{TP} + \operatorname{FP}, \operatorname{TP} + \operatorname{FN}\) と異なっているので，調和平均を用いて以下のように計算します：
\[ \begin{align*} \operatorname{F1-score} &= \frac{1}{\frac{1}{2} \left(\frac{1}{\text{precision}} + \frac{1}{\text{recall}}\right)}\\ &= \frac{2}{\frac{1}{\text{precision}} + \frac{1}{\text{recall}}} \end{align*} \]
なお，これはウェイトが等しい場合を意味しており，weighted harmonic meanへ拡張する場合は以下のように \(\operatorname{F_\beta-score}\) を用いて計算します
\[ \operatorname{F_\beta-score} = \frac{1 + \beta^2}{\frac{1}{\text{precision}} + \frac{\beta^2}{\text{recall}}} \]
ウェイトが \(\frac{1}{1 + \beta^2}, \frac{\beta^2}{1 + \beta^2}\) の形を取っているのは一見不自然に見えますが，その考察で面白いのがvan Rijs-bergen’s E (effectiveness) functionに基づいた説明です．
Theorem 1.1
定義域が \(\mathbb R_+\) の確率変数 \(X\) を考える（つまり \(X > 0\)）. このとき,
として以下が常に成り立つ
\[ H_x \leq G_x \leq \overline{X} \]
Def: median
データ \(x_1, \cdots, x_n\) を小さい順に並び替えた順序統計量
\[ x_{[1]} < \cdots < x_{[n]} \]
について，真ん中の値を中央値という．つまり，
\[ \operatorname{Med}(X) = \left\{\begin{array}{cl} x_{[k]} & \text{where } n = 2k-1\\ \displaystyle \frac{x_{[k]} + x_{[k+1]}}{2} & \text{where } n = 2k \end{array}\right. \]
▶ Hodges-Lehmann推定量
標本の中からペアを選び，そのヘアの平均の中央値を用いて中央値を推定するのがホッジスレーマン推定値です．
\[ \hat\mu_{HL} = \operatorname{Med}\bigg(\bigg\{\frac{x_i + x_j}{2}\bigg\}_{1\leq i \leq j \le n}\bigg) \]
computation上少し重たいですが計算例として以下，
import itertools
def HL_mean(x: list[tuple]):
return np.median([np.mean(t) for t in itertools.combinations(x, 2)])
X_0 = np.array([5.6, 5.7, 5.4, 5.5, 5.8, 5.2, 5.3, 5.6, 5.4, 55.5])
print(HL_mean(X_0))5.55Theorem 1.2 : L1距離最小化推定量としての中央値
標本 \(X_1, \cdots, X_n\) についてmedianを以下のように定義する
\[ \operatorname{Med}(X) = \left\{\begin{array}{cl} x_{[k]} & \text{where } n = 2k-1\\ \displaystyle \frac{x_{[k]} + x_{[k+1]}}{2} & \text{where } n = 2k \end{array}\right. \]
このとき，medianは次の性質を満たす
\[ \operatorname{Med}(X) = \arg\min_a \sum \vert X_i - a\vert \]
Theorem 1.3 : L1距離最小化推定量としての中央値 - 連続確率変数の場合
\(F(x)\) を分布関数に持つ連続確率変数 \(X\) について，
\[ \operatorname{Med}(X) = \arg\min_a \mathbb E[\vert X - a\vert] \]
が成り立つ．
Theorem: Asymptotic distribution of sample quantile-p
\(y_1, \cdots, y_n\) を density function \(f\) 及びquantile function \(Q^{(p)}\) を持つ分布からのi.i.dとします．このとき, sample quantile \(\hat Q^{(p)}\) は
\[ \sqrt{n}(\hat Q^{(p)} - Q^{(p)}) \rightarrow_d \mathbb N\left(0,\frac{p(1-p)}{f(Q^{(p)})^2}\right) \]
ここでは，連続変数分布を想定して解説します．連続なdensity function \(f_x\) を持つ連続確率分布 \(F\) という分布について
\[ X = \{x_1, \cdots, x_n\} \overset{\mathrm{iid}}{\sim} F \]
と確率変数列が与えられたとします．この確率変数に対して
\[ Z_i\equiv 1\{x_i \leq x\} \]
という変数を考えます．この \(Z_i\) は Bernoulli分布に従うと考えられるので，\(F\) のCDFを \(F_X\) とおくと
\[ \begin{align*} \mathbb E(Z_i) &= \mathbb E\left(I\{X_i\le x\}\right) = P(X_i\le x)=F_X(x)\\ \operatorname{Var}(Z_i) &= F_X(x)[1-F_X(x)] \end{align*} \]
ここで，\(Z_i\) のsample meanを以下のように定義する．
\[ Y_n(x) = \frac 1n\sum_{i=1}^nZ_i \]
このように定義した \(Y_n(x)\) はいわゆる経験分布関数 \(F_n(x)\) であるとみなせます．また，定義より
\[ \begin{align*} &E[F_n(x)] = F_X(x)\\ &\operatorname{Var}(F_n(x)) = (1/n)F_X(x)[1-F_X(x)]\\ &\sqrt n\Big(F_n(x) - F_X(x)\Big) \rightarrow_d \mathbb N\left(0,F_X(x)[1-F_X(x)]\right) \because{\text{CLT}} \end{align*} \]
ここでCDFの逆関数 \(F^{-1}_X\) とする(monotonicityより自明)と delta methodを用いると
\[ \begin{align*} &\frac {d}{dt}F^{-1}_X(t) = \frac 1{f_x\left(F^{-1}_X(t)\right)}\\ &\sqrt n\Big(F^{-1}_X(F_n(x)) - F^{-1}_X(F_X(x))\Big) \rightarrow_d \mathbb N\left(0,\frac {F_X(x)[1-F_X(x)]}{\left[f_x\left(F^{-1}_X(F_X(x))\right)\right]^2} \right) \end{align*} \]
つまり，
\[ \sqrt n\Big(F^{-1}_X(F_n(x)) - x\Big) \rightarrow_d \mathbb N\left(0,\frac {F_X(x)[1-F_X(x)]}{\left[f_x(x)\right]^2} \right) \]
ここで \(x = m\)(population median)と設定すると
\[ \sqrt n\Big(F^{-1}_X(F_n(m)) - m\Big) \rightarrow_d \mathbb N\left(0,\frac {1}{\left[2f_x(m)\right]^2} \right) \]
また，
\[ F^{-1}_X(\hat F_n(m)) = \inf\{x : F_X(x) \geq \hat F_n(m)\} = \inf\{x : F_X(x) \geq \frac 1n \sum_{i=1}^n I\{X_i\leq m\}\} \]
より, 不等式のRHSは 1/2 に収束するので \(F^{-1}_X(\hat F_n(m))\) はsample mean \(\hat m\) に収束することがわかる．従って，
\[ \sqrt n\Big(\hat m - m\Big) \rightarrow_d \mathbb N\left(0,\frac {1}{\left[2f_x(m)\right]^2} \right) \]
Theorem 1.4 標本平均とMedianの距離
独立に同一の分布に従う確率変数列 \(\{X_i\}_{i=1}^n\) を考える．この標本平均を \(\overline{X}\), 不偏分散を \(S^2\), メディアン \(X_{\operatorname{med}}\) とすると
\[ \vert \overline{X} - X_{\operatorname{med}}\vert < S \]
Def: Variance
mean \(\mu\) をもつ確率変数 \(X\) の分散は
\[ \begin{align*} \operatorname{Var}(X) &= \mathbb E[(X - \mu)^2]\\ &= \mathbb E[X^2] - \mathbb E[X]^2 \end{align*} \]
標準偏差(standard deviation)は分散のsquare rootで定義される．
上の定義より標準偏差について以下のことがわかります：
Theorem 1.5
確率変数 \(X\sim D(\mu, \sigma^2)\) とする．このとき
\[ \begin{align*} \mathbb E[\vert X - \mu\vert] \leq \sigma \end{align*} \]
確率変数のバラツキを表す指標として範囲(Range), 四分位範囲（IQR）, 分散がありますが，分布の中心の位置が著しく異なるような場合 には，これらを用いて分布の散らばり具合を比較することは難しいです．このような場合に，変動係数(Coefficient of variation) という単位のない統計量を用いたりします．
Def: Coefficient of variation
比例尺度にもとづく測定が行われ，その測定結果をnon-negative 確率変数 \(X_i >0\) で表すとする．\(\{X_i\}_{i=1}^n\) の標本平均を \(\overline{X}\), 標本標準偏差を \(S\) としたとき，変動係数は以下のように計算される
\[ \operatorname{C_V} = \frac{S}{\overline{X}} \]
上記の定義より以下のことがわかります
Example 1.1 東京 日平均気温の月平均値（℃）
国土交通省気象庁より以下のように東京都の気象データを取得します．
import polars as pl
import plotly.express as px
from regmonkey_style import stylewizard as sw
sw.set_templates("regmonkey_line")
data = {
"Year": [2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023],
"Jan": [5.8, 6.1, 5.8, 4.7, 5.6, 7.1, 5.4, 4.9, 5.7],
"Feb": [5.7, 7.2, 6.9, 5.4, 7.2, 8.3, 8.5, 5.2, 7.3],
"Mar": [10.3, 10.1, 8.5, 11.5, 10.6, 10.7, 12.8, 10.9, 12.9],
"Apr": [14.5, 15.4, 14.7, 17.0, 13.6, 12.8, 15.1, 15.3, 16.3],
"May": [21.1, 20.2, 20.0, 19.8, 20.0, 19.5, 19.6, 18.8, 19.0],
"Jun": [22.1, 22.4, 22.0, 22.4, 21.8, 23.2, 22.7, 23.0, 23.2],
"Jul": [26.2, 25.4, 27.3, 28.3, 24.1, 24.3, 25.9, 27.4, 28.7],
"Aug": [26.7, 27.1, 26.4, 28.1, 28.4, 29.1, 27.4, 27.5, 29.2],
"Sep": [22.6, 24.4, 22.8, 22.9, 25.1, 24.2, 22.3, 24.4, 26.7],
"Oct": [18.4, 18.7, 16.8, 19.1, 19.4, 17.5, 18.2, 17.2, 18.9],
"Nov": [13.9, 11.4, 11.9, 14.0, 13.1, 14.0, 13.7, 14.5, 14.4],
"Dec": [9.3, 8.9, 6.6, 8.3, 8.5, 7.7, 7.9, 7.5, 9.4],
}
df = pl.DataFrame(data)
df_unpivoted = df.unpivot(index="Year", variable_name="Month", value_name="Celsius")
fig = px.line(df_unpivoted, x="Month", y="Celsius", color="Year", title="東京都月別平均気温")
fig.show()ここで，以下のようにFahrenheitに変換してみます．
df_unpivoted = df_unpivoted.with_columns(
(pl.col("Celsius") * 9 / 5 + 32).alias("Fahrenheit")
)
df_unpivoted.head()| Year | Month | Celsius | Fahrenheit |
|---|---|---|---|
| i64 | str | f64 | f64 |
| 2015 | "Jan" | 5.8 | 42.44 |
| 2016 | "Jan" | 6.1 | 42.98 |
| 2017 | "Jan" | 5.8 | 42.44 |
| 2018 | "Jan" | 4.7 | 40.46 |
| 2019 | "Jan" | 5.6 | 42.08 |
ここで，Celsius, Fahrenheit両方のカラムについて変動係数を計算します．
def compute_cv(df, col: str) -> float:
return df_unpivoted[col].std() / df_unpivoted[col].mean()
print(
"""Celsius CV: {:.2f}, Fahrenheit CV: {:.2f}""".format(
compute_cv(df_unpivoted, "Celsius"), compute_cv(df_unpivoted, "Fahrenheit")
))Celsius CV: 0.45, Fahrenheit CV: 0.22このように，Celsius, Fahrenheitともに同じ温度を単位の違う方法で表しているのにも関わらず変動係数は異なります．Celsius, Fahrenheitともに間隔尺度であること，及び変数変換の観点からもlocation/scale parameterを異なる値で調整しているので同じデータを扱っているにも関わらずCVが一致しないという現象が発生してしまいます．
▶ Bias Correction
サイズ \(N\) のsampleベースで計算された変動係数はpopulation変動係数 \(\gamma_V\) と比較して過小推定されているということが知られています．population変動係数のunbiased estimate \(\widehat{\operatorname{C_V}}\) は以下のように計算されます
\[ \widehat{\operatorname{C_V}} = \left(1 + \frac{1}{4N}\right)\operatorname{C_V} \]
▶ \(\sqrt{n}\) の法則
確率変数 \(X_1, \cdots, X_n\) 独立に同一の分布に従うとし，これらの期待値と分散を \(\mu, \sigma^2\) とします．
このとき，\(\tilde X_n = \sum^n_{i=1} X_i\) とすると，
\[ \begin{align*} \mathbb E[\tilde X_n] &= n\mu\\ \operatorname{Var}(\tilde X_n) &= n\sigma^2\\ \operatorname{std}(\tilde X_n) &= \sqrt{n}\sigma \end{align*} \]
期待値は \(n\) のオーダーで増える一方，標準偏差は \(\sqrt{n}\) オーダーなので
\[ \begin{align*} \operatorname{CV}(\tilde X_n) = \frac{\sqrt{n}\sigma}{n\mu} = \frac{\sigma}{\sqrt{n}\mu} \end{align*} \]
ここから，平均と標準偏差の比率は \(n\) が増えるほど，小さくなる = データのバラツキが相対的に小さくなることがわかります．
とある店AとBについて，それぞれの各月の来店者人数は各店舗が立地している地域の人口に比例すると解釈して， \(X_A\sim\operatorname{Binom}(200, 0.5), Y_B\sim\operatorname{Binom}(1000, 0.5)\) である場合を考えます．
それぞれの店舗について来月来店人数を予測したいとき，\(p_A = p_B = 0.5, N_A= 200, N_B= 1000\) であると判断し，その下で期待値をとって \(100, 500\) と予測したとします． このとき，予測値を基準に予測誤差 5% を上回りそうな確率はそれぞれ
\[ \begin{align*} 1 - \Pr(\vert X_A - 100\vert < 5) &\approx 0.481\\ 1 - \Pr(\vert X_B - 500\vert < 25) &\approx 0.114 \end{align*} \]
と店舗 B のほうがMAPE的により良い予測精度が得られそうなことがわかります．母集団parameterもそれぞれ正しく判断できているにも関わらず, MAPEという指標では店舗 B のほうが相対的に良さそうな結果がでてしまうと推察できます．