チェビシェフの不等式の導出とシミュレーション

チェビシェフの不等式の導出とシミュレーション

チェビシェフの不等式の導出メモ．

チェビシェフの不等式とは

\begin{eqnarray} P(|X - \mu| \geq a \sigma) \leq \frac{1}{a^2} \end{eqnarray} で表される．

確率変数 \(X\) が平均から標準偏差の\(a\)倍離れている確率は \( \frac{1}{a^2} \) 以下であることを表している．

例えば，\( a=2 \) の場合， 平均から標準偏差の\(2\)倍離れている確率は \( \frac{1}{2^2} \) 以下である．

別の書き方をすれば， \begin{eqnarray} 1 - \frac{1}{a^2} \leq P(|X - \mu| \leq a\sigma) \end{eqnarray} で表される．

この場合， 確率変数 \(X\) が平均から標準偏差の\(a\)倍以内の範囲にある確率は \( \frac{1}{a^2} \) 以上であることを表している．

例えば，\( a=2 \) の場合， 平均から標準偏差の\(2\)倍以内の確率は \( 1 - \frac{1}{2^2} = \frac{3}{4} \) 以上になる．

利点

不等式ということで，ざっくりしかわからないが， チェビシェフの不等式は全ての分布について成り立つ．

なので，平均と標準偏差さえわかれば， 謎の分布に従っていても，平均から標準偏差の2倍離れているデータが得られたとき，確率\( \frac{1}{4} \) 以下の珍しいことが起こったと概算できる．

チェビシェフの不等式の導出

連続変数の場合の分散の定義式で積分を3つの範囲に分ける． ただし，\(a > 0, f(x) \geq 0\) とする． \begin{eqnarray*} \sigma^2 &=& \int_{-\infty}^{\infty} (x-\mu^2) f(x) dx \\ &=& \int_{-\infty}^{\mu - a\sigma} (x-\mu^2) f(x) dx +\int_{\mu - a\sigma}^{\mu + a\sigma} (x-\mu^2) f(x) dx +\int_{\mu + a\sigma}^{\infty} (x-\mu^2) f(x) dx \end{eqnarray*}

右辺第2項を \( 0 \) で置き換えると，不等式になり． \begin{eqnarray*} \sigma^2 &\geq& \int_{-\infty}^{\mu - a\sigma} (x-\mu^2) f(x) dx +\int_{\mu + a\sigma}^{\infty} (x-\mu^2) f(x) dx \end{eqnarray*} となる．

確率変数 \( X \) の上側，下側範囲を考えると． $$ X \leq \mu - a \sigma \to X - \mu \leq - a \sigma \\ X \geq \mu + a \sigma \to X - \mu \geq a \sigma $$ これをまとめると，\( (X - \mu)^2 \geq a^2 \sigma^2 \) であるから， \( a^2 \sigma^2 \) というより小さい値に置き換えると， \begin{eqnarray*} \sigma^2 &\geq& \int_{-\infty}^{\mu - a\sigma} a^2\sigma^2 f(x) dx +\int_{\mu + a\sigma}^{\infty} a^2\sigma^2 f(x) dx \\ \sigma^2 &\geq& a^2\sigma^2 \left( \int_{-\infty}^{\mu - a\sigma} f(x) dx +\int_{\mu + a\sigma}^{\infty} f(x) dx \right) \\ \sigma^2 &\geq& a^2 \sigma^2 \left( P(X \leq \mu - a\sigma) + P(X \geq \mu - a \sigma) \right) \end{eqnarray*} となる．

\( X \leq \mu - a \sigma \) と \( X \geq \mu + a \sigma \) をまとめると，\( |X - \mu| \geq a \sigma \) となるので， \begin{eqnarray*} \sigma^2 &\geq& a^2 \sigma^2 \left( P(|X - \mu| \geq a\sigma) \right) \\ 1 &\geq& a^2 \left( P(|X - \mu| \geq a\sigma) \right) \\ P(|X - \mu| \geq a\sigma) &\leq& \frac{1}{a^2} \\ \end{eqnarray*}

チェビシェフの不等式のシミュレーション

本当に成り立つのか

• 標準正規分布

• 二項分布

• ポアソン分布

を例にシミュレーションしてみる．

使用環境は以下である．

環境

• conda: 4.2.12
  • Python 3.5.1
  • scipy 0.18.1
  • pandas 0.18.1
  • jupyter 1.0.0

あらかじめ，scipy をimportしておく．

import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt

標準正規分布の場合

\( X \sim N(0, 1) \) として， \( P(|X| \geq 2) \) となる確率は， $$ P(|X| \geq 2) = P(X \geq 2) + P(X \leq -2) = 0.0455 \leq \frac{1}{4} $$ \( P(|X| \geq 3) \) となる確率は， $$ P(|X| \geq 3) = P(X \geq 3) + P(X \leq -3) = 0.0027 \leq \frac{1}{9} $$ 以上より，チェビシェフの不等式を満たしている．

実際に乱数を1000個発生させて，標準偏差の2倍3倍を超える割合を計算してみる.

  %matplotlib inline

  # simulation
  # N(0, 1)
  norm2sigma_prob = sp.stats.norm.sf(loc=0, scale=1, x=2*1) + (1 - sp.stats.norm.sf(loc=0, scale=1, x=-2*1))
  print("N(0, 1): 2sigma over prob", norm2sigma_prob, "<= 1/4?:", norm2sigma_prob <= 1/2**2)
  norm3sigma_prob = sp.stats.norm.sf(loc=0, scale=1, x=3*1) + (1 - sp.stats.norm.sf(loc=0, scale=1, x=-3*1))
  print("N(0, 1): 3sigma over prob", norm3sigma_prob, "<= 1/9?:", norm3sigma_prob <= 1/3**2)
  # random
  norm = stats.norm.rvs(size=1000, loc=0, scale=1)
  norm2sigma = norm[sp.absolute(norm)>=2*1] # 標準偏差の2倍以上の値を抽出
  print("Random: 2sigma over ratio", norm2sigma.size/1000, "<= 1/4?:", norm2sigma.size/1000 <= 1/2**2)
  norm3sigma = norm[sp.absolute(norm)>=3*1] # 標準偏差の3倍以上の値を抽出
  print("Random: 3sigma over ratio", norm3sigma.size/1000, "<= 1/9?:", norm3sigma.size/1000 <= 1/3**2)

  # plot
  # N(0, 1)
  x_ = sp.linspace(-4, 4, 100)
  plt.plot(x_, stats.norm.pdf(x_), "-r", label="N(0, 1)")
  # hist
  bins = sp.linspace(-4, 4, 12)
  plt.hist(norm, bins=bins, normed=True, alpha=0.6, label="random freq")
  # mu = 0
  plt.vlines(0, 0, 0.5)
  plt.text(-0.2, -0.03, "mu")
  # fill_betweenで塗りつぶし
  x_ = sp.linspace(2, 4)
  plt.fill_between(x_, 0, stats.norm.pdf(x_), facecolor='r', alpha=0.5)
  x_ = sp.linspace(-4, -2)
  plt.fill_between(x_, 0, stats.norm.pdf(x_), facecolor='r', alpha=0.5)
  plt.text(2, -0.03, "mu+2sigma")
  plt.text(-2-1, -0.03, "mu-2sigma")

  plt.xlim(-4, 4)
  plt.ylim(-0.05, 0.5)
  plt.xlabel("random variable x")
  plt.ylabel("f(x)")
  plt.title("Simulation with Normal Distribution", y=-.25)
  plt.legend()
  plt.grid()
  plt.tight_layout()
  plt.show()

結果は以下のようになり，チェビシェフの不等式を満たしている．

Random: 2sigma over ratio 0.039 <= 1/4?: True
Random: 3sigma over ratio 0.001 <= 1/9?: True

Fig. Plot N(0, 1) and Freq of Random sample(1000)

二項分布の場合

同様にして二項分布の場合， \( X \sim Bin(10, 0.5) \) として， 平均と標準偏差は \( \mu = np = 5\), \( \sigma = \sqrt{np(1-p)} = \sqrt{5 \cdot 0.5} = \sqrt{2.5} \)より， \( P(|X - \mu| \geq 2\sigma) \) となる確率は， $$ P(|X - \mu| \geq 2\sigma) = P(X - \mu \geq 2\sigma) + P(X - \mu \leq -2\sigma) = 0.02148 \leq \frac{1}{4} $$ \( P(|X - \mu| \geq 3\sigma) \) となる確率は， $$ P(|X - \mu| \geq 3 \sigma) = P(X - \mu \geq 3\sigma) + P(X - \mu \leq -3\sigma) =0.00195 \leq \frac{1}{9} $$ 以上より，チェビシェフの不等式を満たしている．

二項分布でも同様にして，標準偏差の2倍3倍を超える割合を計算してみる.

%matplotlib inline

# simulation
# Bin(10, 0.5)
n = 10
p = 0.5
mu = n * p
sigma = sp.sqrt(n * p * (1 - p))

bin2sigma_prob = sp.stats.binom.sf(k=mu+2*sigma, n=n, p=p) + sp.stats.binom.cdf(k=mu-2*sigma, n=n, p=p)
print("Bin(10, 0.5): 2sigma over prob", bin2sigma_prob, "<= 1/4?:", bin2sigma_prob <= 1/2**2)
bin3sigma_prob = sp.stats.binom.sf(k=mu+3*sigma, n=n, p=p) + sp.stats.binom.cdf(k=mu-3*sigma, n=n, p=p)
print("Bin(10, 0.5): 3sigma over prob", bin3sigma_prob, "<= 1/9?:", bin3sigma_prob <= 1/3**2)
# random
bino = stats.binom.rvs(size=1000, n=n, p=p)
bino2sigma = bino[sp.absolute(bino - mu) >= (2 * sigma)] # 標準偏差の2倍以上の値を抽出
print("2sigma over ratio", bino2sigma.size/1000, "<= 1/4?:", bino2sigma.size/1000 <= 1/2**2)
bino3sigma = bino[sp.absolute(bino - mu)>= (3 * sigma)] # 標準偏差の3倍以上の値を抽出
print("3sigma over ratio", bino3sigma.size/1000, "<= 1/9?:", bino3sigma.size/1000 <= 1/3**2)


# plot
# Bin(n, p)
x_ = sp.arange(0, n+1, 1)
plt.plot(x_, stats.binom.pmf(x_, n=n, p=p), "-r", label="Bin(10, 1/2)")
# hist
bins = sp.linspace(0, n, n)
plt.hist(bino, bins=bins, normed=True, alpha=0.6, label="random freq")
# mu = 0
plt.vlines(mu, 0, 0.3)
plt.text(mu-0.2, -0.02, "mu")
# fill_betweenで塗りつぶし
x_ = sp.arange(round(mu + 2*sigma), n+1, 1)
plt.fill_between(x_, 0, stats.binom.pmf(k=x_, n=n, p=p), facecolor='r', alpha=0.5)
plt.text(mu+2*sigma, -0.03, "mu+2sigma")
x_ = sp.arange(0, round(mu - 2*sigma)+1, 1)
plt.fill_between(x_, 0, stats.binom.pmf(k=x_, n=n, p=p), facecolor='r', alpha=0.5)
plt.text(mu-2*sigma-1, -0.03, "mu-2sigma")

plt.xlabel("random variable x")
plt.ylabel("f(x)")
plt.title("Simulation with Binom Distribution", y=-.25)
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()

結果は以下のようになり，チェビシェフの不等式を満たしている．

  2sigma over ratio 0.024 <= 1/4?: True
  3sigma over ratio 0.001 <= 1/9?: True

Fig. Plot Bin(10, 0.5) and Freq of Random sample(1000)

Poisson分布の場合

同様にしてPoisson分布の場合， \( X \sim Po(10) \) として， \( P(|X - \mu| \geq 2\sigma) \) となる確率は， $$ P(|X - \mu| \geq 2\sigma) = P(X - \mu \geq 2\sigma) + P(X - \mu \leq -2\sigma) = 0.0374 \leq \frac{1}{4} $$ \( P(|X - \mu| \geq 3\sigma) \) となる確率は， $$ P(|X - \mu| \geq 3 \sigma) = P(X - \mu \geq 3\sigma) + P(X - \mu \leq -3\sigma) = 0.0035 \leq \frac{1}{9} $$ 以上より，チェビシェフの不等式を満たしている．

Poisson分布でも同様にして，計算する．

%matplotlib inline

# simulation
# Po(10)
mu = 10
sigma = sp.sqrt(mu)

po2sigma_prob = stats.poisson.sf(k=mu+2*sigma, mu=mu) + stats.poisson.cdf(k=mu-2*sigma, mu=mu)
print("Po(10): 2sigma over prob", po2sigma_prob, "<= 1/4?:", po2sigma_prob <= 1/2**2)
po3sigma_prob = stats.poisson.sf(k=mu+3*sigma, mu=mu) + stats.poisson.cdf(k=mu-3*sigma, mu=mu)
print("Po(10): 3sigma over prob", po3sigma_prob, "<= 1/9?:", po3sigma_prob <= 1/3**2)
# random
po = stats.poisson.rvs(size=1000, mu=mu)
po2sigma = po[sp.absolute(po - mu) >= (2 * sigma)] # 標準偏差の2倍以上の値を抽出
print("Random: 2sigma over ratio", po2sigma.size/1000, "<= 1/4?:", po2sigma.size/1000 <= 1/2**2)
po3sigma = po[sp.absolute(po - mu)>= (3 * sigma)] # 標準偏差の3倍以上の値を抽出
print("Random: 3sigma over ratio", po3sigma.size/1000, "<= 1/9?:", po3sigma.size/1000 <= 1/3**2)


# plot
# Po
x_ = sp.arange(0, n*2, 1)
plt.plot(x_, stats.poisson.pmf(x_, mu=mu), "-r", label="Po(10)")
# hist
bins = sp.linspace(0, n*2, n*2)
plt.hist(po, bins=bins, normed=True, alpha=0.6, label="random freq")
# mu = 10
plt.vlines(mu, 0, 0.2)
plt.text(mu-0.3, -0.02, "mu")
# fill_betweenで塗りつぶし
x_ = sp.arange(round(mu + 2*sigma), n*2+1, 1)
plt.fill_between(x_, 0, stats.poisson.pmf(k=x_, mu=mu), facecolor='r', alpha=0.5)
plt.text(mu+2*sigma, -0.02, "mu+2sigma")
x_ = sp.arange(0, round(mu - 2*sigma)+1, 1)
plt.fill_between(x_, 0, stats.poisson.pmf(k=x_, mu=mu), facecolor='r', alpha=0.5)
plt.text(mu-2*sigma-2, -0.02, "mu-2sigma")

plt.ylim(0, 0.2)
plt.xlabel("random variable x")
plt.ylabel("f(x)")
plt.title("Simulation with Poisson Distribution", y=-.25)
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()

結果は以下のようになり，チェビシェフの不等式を満たしている．

Random: 2sigma over ratio 0.037 <= 1/4?: True
Random: 3sigma over ratio 0.002 <= 1/9?: True

Fig. Plot Po(10) and Freq of Random sample(1000)

参考文献

• 薩摩順吉, "確率・統計", 岩波書店, p.46

  確率・統計 (理工系の数学入門コース 7)