確率変数の不等式を視覚的に理解する〜マルコフの不等式・チェビシェフの不等式〜

  • Wikipediaのこのページに沿って\mu(\{x \in X ; f(x) \ge t\}) \le \frac{1}{g(t)}\int_{f(x) \ge t} g(f(x)) d\muとはどういうことかを順を追って視覚的に確認してみます

  • 以下はRmdファイルです
---
title: "確率変数の不等式を視覚的に理解する〜マルコフの不等式・チェビシェフの不等式〜"
author: "ryamada"
date: "Thursday, May 14, 2015"
output: html_document
---


# はじめに

確率分布の平均値や分散がわかると、その確率変数がある範囲の値を取る確率がどのくらいになるかについてわかることがある。

このわかり方として、「これこれの範囲の値を取る確率はこれこれ以下である」というように、不等式を用いた形でわかることがあり、そのような不等式を、確率論の不等式と言い、いくつも知られている。

ビッグデータ活用にあたって、全部のデータを使うわけには行かないが、この範囲に収まっているはずだ、ということを明言する場合にも使われる。

そんな確率変数の不等式のうち、ごく基本的な二つ、マルコフの不等式とチェビシェフの不等式を、図を使って視覚的に理解することを目指す。


# ベータ分布を例に

ベータ分布を描いてみる。

```{r}
a <- 3
b <- 2
x <- seq(from=0,to=1,length=1000)
y <- dbeta(x,a,b)
plot(x,y,type="l")
```

このように$x \in X=[0,1]$でその確率密度は$y \ge 0$の分布である。

今、$x \ge 0.8$である確率を$\mu(\{x \in X; x \ge 0.8\})$ と表すことにすれば、これは

```{r}
plot(x,y,type="l")
abline(v=0.8,col=2)
p0.8 <- pbeta(0.8,a,b)
1-p0.8
```
で計算できる。ベータ分布の累積分布関数の値を1から引いた値になっている。

累積分布の値を0-1について密に調べれば

```{r}
t <- seq(from=0,to=1,length=1000)
z <- pbeta(t,a,b)
plot(t,z,type="l")
abline(v=0.8,col=2)
abline(h=p0.8,col=3)
```

のようになる。

興味のある$\mu(\{x \in X; x \ge 0.8\})$というのは、累積分布
関数について、閾値0.8に縦線()を引き、曲線との交点を通る横線()を引いて、その横線の高さを1から引いた長さのことである。

この累積分布グラフの縦軸と横軸を入れ替えると次のようになる。

```{r}
plot(z,t,type="l")
abline(h=0.8,col=2)
abline(v=p0.8,col=3)
```


このグラフで言えば、曲線が閾値0.8以上になっている範囲の$z$の軸の長さが知りたい値である。

この『長さを測って』もよいけれど、少し視点を変えて、次のような関数を作って、『積分』して長さと同じ数値を得ることにする。

```{r}
w <- rep(0,length(z))
w[which(t>=0.8)] = 1
matplot(z,cbind(t,w),type="l",col=c(1,4))
abline(h=0.8,col=2)
abline(v=p0.8,col=3)
```

付け加えた青い折れ線は、$t < 0.8$では0をとり$t \ge 0.8$では1を取る関数のそれである。

このとき、$mu$はこの青い折れ線グラフの積分値に等しい(水色の長方形)。

```{r}
w <- rep(0,length(z))
w[which(t>=0.8)] = 1
matplot(z,cbind(t,w),type="l",col=c(1,4))
abline(h=0.8,col=2)
abline(v=p0.8,col=3)
rect(p0.8,0,1,1,col=5)
```

ここで、紫の長方形を考える。この面積は、水色の面積($\mu$)0.8倍であるから、紫の長方形の面積がわかれば、$\mu$の面積は計算できる。


```{r}
w <- rep(0,length(z))
w[which(t>=0.8)] = 1
matplot(z,cbind(t,w),type="l",col=c(1,4))
abline(h=0.8,col=2)
abline(v=p0.8,col=3)
rect(p0.8,0,1,1,col=5)
rect(p0.8,0,1,0.8,col=6)
```

水色の長方形を取り除いて表示しなおしてみる。

```{r}
w <- rep(0,length(z))
w[which(t>=0.8)] = 1
plot(z,t,type="h")
#matplot(z,cbind(t,w),type="h",col=c(1,4))
abline(h=0.8,col=2)
abline(v=p0.8,col=3)
#rect(p0.8,0,1,1,col=5)
rect(p0.8,0,1,0.8,col=6)
```

この紫の長方形の面積は、黒い曲線の積分値以下であることは見て取れる。
したがって、黒い曲線の積分値がわかればよい。
この面積は、もう一度縦軸と横軸を入れ替えなおしてやって考えれば
```{r}
plot(t,z,type="h")
```
のようになる。

この下面積を計算するにあたり、縦軸は確率変数の値、横軸の増分は確率密度に相当する値なので

$$
\int_0^1 x dPr
$$
となる。これは、確率変数の値を生起確率で重みづけした平均の算出式になっており、このベータ分布の期待値Eになる。

結局、興味を持っている$\mu$は閾値0.8とEとを使って

$$
0.8 \mu \le E
$$
と表せることがわかり
$$
\mu \le \frac{E}{0.8}
$$
である。

実際、ベータ分布のEは計算できて
```{r}
E <- a/(a+b)
E
```
確かに、理論的な$\mu$値と上記の方法で求めた値との関係は成り立っている。
```{r}
1-p0.8 < E/0.8
```

これを色々な閾値について実施すれば以下のようになる。
横軸が閾値で、
黒線が理論値、赤線が計算した上限値である。
右のパネルはy軸のゼロ付近を拡大表示したものである。

```{r}
up.bd <- E/t
par(mfcol=c(1,2))
matplot(t,cbind(1-z,up.bd),type="l")
matplot(t,cbind(1-z,up.bd),type="l",ylim=c(0,3))
par(mfcol=c(1,1))
```

# Markovの不等式

もう一度、紫の長方形と黒い領域の図を見直してみる。
```{r}
w <- rep(0,length(z))
w[which(t>=0.8)] = 1
plot(z,t,type="h")
#matplot(z,cbind(t,w),type="h",col=c(1,4))
abline(h=0.8,col=2)
abline(v=p0.8,col=3)
#rect(p0.8,0,1,1,col=5)
rect(p0.8,0,1,0.8,col=6)
```

紫の面積が黒を積分した面積以下になることが自明であるためには、黒領域を決める曲線が正の値を取っている必要がある。
なぜなら、曲線が負の値もとると、その部分の積分によって「面積」は小さくなってしまうからである。

したがって、ベータ分布のように確率変数が取る値が0以上であることを条件とすれば、確率分布の期待値を用いて、その確率変数が閾値以上を取る確率の上限を与えるのが、Markovの不等式である。

$$
\mu(\{x \in X; f(x) \ge t\}) \le \frac{E}{t}
$$

# 負の値もとる確率変数ではうまく行かない

負の値もとる確率変数ではMarkovの不等式がうまくいかないことを示して億。

今度はベータ分布ではなく正規分布を用いる。
ベータ分布では確率変数は0-1の値を取ったが正規分布では実数全体を取る。

以下の例では、正規分布の期待値が負になるようにパラメタaを負にしている。

分布を描く。
```{r}
a <- -3
b <- 2
x <- seq(from=-20,to=20,length=1000)
y <- dnorm(x,a,b)
plot(x,y,type="l")
```

閾値を描きいれ、理論的な$\mu$値を計算する。
```{r}
thres <- 0.8
plot(x,y,type="l")
abline(v=thres,col=2)
p0.8 <- pnorm(thres,a,b)
1-p0.8
```
累積分布を描く。

```{r}
t <- seq(from=-20,to=20,length=1000)
z <- pnorm(t,a,b)
plot(t,z,type="l",ylim=c(0,1))
abline(v=thres,col=2)
abline(h=p0.8,col=3)
```

この累積分布グラフの縦軸と横軸を入れ替えると次のようになる。

横軸はベータ分布のときと同じで0-1の範囲であるが、縦軸は、実数全体になる。

```{r}
plot(z,t,type="l")
abline(h=thres,col=2)
abline(v=p0.8,col=3)
```

階段状の関数が同じように描ける。

```{r}
w <- rep(0,length(z))
w[which(t>=thres)] = 1
matplot(z,cbind(t,w),type="l",col=c(1,4))
abline(h=thres,col=2)
abline(v=p0.8,col=3)
```

長方形も同様に描ける。

```{r}
w <- rep(0,length(z))
w[which(t>=thres)] = 1
matplot(z,cbind(t,w),type="l",col=c(1,4))
abline(h=thres,col=2)
abline(v=p0.8,col=3)
rect(p0.8,0,1,1,col=5)
```


```{r}
w <- rep(0,length(z))
w[which(t>=thres)] = 1
matplot(z,cbind(t,w),type="l",col=c(1,4))
abline(h=thres,col=2)
abline(v=p0.8,col=3)
rect(p0.8,0,1,1,col=5)
rect(p0.8,0,1,thres,col=6)
```

```{r}
w <- rep(0,length(z))
w[which(t>=thres)] = 1
plot(z,t,type="h")
#matplot(z,cbind(t,w),type="h",col=c(1,4))
abline(h=thres,col=2)
abline(v=p0.8,col=3)
#rect(p0.8,0,1,1,col=5)
rect(p0.8,0,1,thres,col=6)
```


```{r}
plot(t,z,type="h")
```

実際、正規分布のEはパラメタで指定されていて
```{r}
E <- a
E
```
理論的な$\mu$値と上記の方法で求めた値との関係が成り立っていない。
```{r}
1-p0.8 < E/thres
```
次のグラフを見ると、E(積分値)が負になるため、閾値が正の範囲で「上限値」のはずの赤線が負になっており、上限としてふさわしくないことが見て取れる。
```{r}
up.bd <- E/t
par(mfcol=c(1,2))
matplot(t,cbind(1-z,up.bd),type="l")
matplot(t,cbind(1-z,up.bd),type="l",ylim=c(-1.1,1.1))
par(mfcol=c(1,1))
```

# 正の値を取らせるように拡張する

正規分布では負の値を取るのでMarkov不等式が使えなかった。
正規分布のような負の値を取る分布において、「正の値」だけを取るように考えることができるだろうか。

一つは、確率変数が正の値を取る範囲だけを対象にして考えるという手がある。

もう一つは、絶対値を取って負の値も正の値にしてしまうという手である。
例えば、平均からのずれを問題にし、ずれの方向は気にしないというのであれば、そのようにすることができる。

では、前の例の正規分布で、平均からのずれの関数で同じことをやってみる。

```{r}
a <- -3
b <- 2
x <- seq(from=-20,to=20,length=1000)
y <- dnorm(x,a,b)
plot(x,y,type="l")
```

閾値を描きいれ、理論的な$\mu$値を計算する。
```{r}
thres <- 0.8
#plot(x,y,type="l")
#abline(v=thres,col=2)
p0.8 <- pnorm(c(-thres,thres),0,b)
#1-p0.8
```
累積分布を描く。

```{r}
t <- seq(from=-20,to=20,length=1000)
z <- pnorm(t,a,b)
plot(t,z,type="l",ylim=c(0,1))
#abline(v=thres,col=2)
#abline(h=p0.8,col=3)
```

この累積分布グラフの縦軸と横軸を入れ替えると次のようになる。

横軸はベータ分布のときと同じで0-1の範囲であるが、縦軸は、実数全体になる。

```{r}
plot(z,t,type="l")
#abline(h=thres,col=2)
#abline(v=p0.8,col=3)
```

ここで、平均からのずれを問題にするので、縦軸の値は、平均aの場合に0になるようにシフトし、そのうえで、縦軸の値が負の部分を正の側に折り返すと、以下のように、平均からの距離が閾値0.8である範囲が図示できる。

```{r}
abs.d <- abs(t-a)
matplot(z,cbind(t,abs.d),type="l")
abline(h=thres,col=2)
abline(v=p0.8,col=3)
```
この平均からの距離に対する
階段状の関数が同じように描ける。

```{r}
w <- rep(0,length(z))
w[which(abs.d>=thres)] = 1
matplot(z,cbind(abs.d,w),type="l",col=c(1,4))
abline(h=thres,col=2)
abline(v=p0.8,col=3)
```

長方形も同様に描ける。

```{r}
w[which(abs.d>=thres)] = 1
matplot(z,cbind(abs.d,w),type="l",col=c(1,4))
abline(h=thres,col=2)
abline(v=p0.8,col=3)
rect(p0.8[1],0,0,1,col=5)
rect(p0.8[2],0,1,1,col=5)
```


```{r}
w[which(abs.d>=thres)] = 1
matplot(z,cbind(abs.d,w),type="l",col=c(1,4))
abline(h=thres,col=2)
abline(v=p0.8,col=3)
rect(p0.8[1],0,0,1,col=5)
rect(p0.8[2],0,1,1,col=5)
rect(p0.8[1],0,0,thres,col=6)
rect(p0.8[2],0,1,thres,col=6)
```

```{r}
w <- rep(0,length(z))
w[which(t>=thres)] = 1
plot(z,abs.d,type="h")
#matplot(z,cbind(t,w),type="h",col=c(1,4))
abline(h=thres,col=2)
abline(v=p0.8,col=3)
#rect(p0.8,0,1,1,col=5)
rect(p0.8[1],0,0,thres,col=6)
rect(p0.8[2],0,1,thres,col=6)
```

この図を見ると、紫の長方形の面積は、黒の積分面積以下であることがわかり、正規分布の場合であっても、平均からの距離が閾値以上である確率の上限が計算できることがわかる。

しかしながら、この黒の領域の面積は、Markovの不等式の場合に期待値Eであったような分布の簡単な期待値ではない。

# Chebyshevの不等式

ここで、黒領域を決める曲線を$|f(x)-E|$から、$(|f(x)-E|^2)$に変更してやる。

図にすれば以下の通りである。

閾値も二乗($0.8^2$)ので、その線も加えておく。
```{r}
sq.d <- abs.d^2
matplot(z,cbind(abs.d,sq.d),type="l",col=c(1,6),ylim=c(0,10))
abline(h=thres,col=2)
abline(h=thres^2,col=2)
abline(v=p0.8,col=3)
```

平均からの距離の関数を平均からの距離の二乗の関数に取り換えても、上限を与えることはできそうです。
そして生起確率で重みづけした平均からのずれの二乗の和の平均は分散のことなので、分布の分散がわかれば、上限が計算できる。




実際、正規分布のEはパラメタで指定されていて
```{r}
V <- b^2
V
```
理論的な$\mu$値と上記の方法で求めた値との関係が成り立っていない。
```{r}
1-p0.8 < V/(thres^2)
```

```{r}
# 平均からの距離の閾値
t2 <- seq(from=0,to=20,length=1000)
up.bd <- V/(t2^2)
# sd=bの正規分布に照らして閾値以上、平均からずれる確率
v.more <- (1-pnorm(t2,0,b))*2
par(mfcol=c(1,2))
matplot(t,cbind(v.more,up.bd),type="l")
matplot(t,cbind(v.more,up.bd),type="l",ylim=c(-1.1,1.1))
par(mfcol=c(1,1))
```

このように
$$
\mu(\{x \in X; |x-E| \ge t\}) \le \frac{V}{t^2}
$$
という不等式関係がある。
これをChebyshevの不等式と言う。

# 一般化する

はじめにMarkovの不等式を見た。変数が非負でないとうまくいかなかった。
次に変数が負であっても絶対値を取るなどして、非負にしてやれば、不等式があることがわかった。その例としてChebyshevの不等式があり、確率変数の分散を使って期待値からのずれが閾値より大きい割合の上限が示せることがわかった。

この2つの不等式の導き方を一般化するなら、何か関数があるときに、それを何かしらの方法で非負に変換してやることで、変換後の関数の値が閾値以上である割合の上限を与えることができる。

今、この非負化変換として、元の関数の値の順序を変えないような変換を選べば、元の関数の値が閾値以上である割合の上限が得られる。

Chebyshevの不等式の場合には、絶対値を取る関数を使ったので、変換の前後で関数の値の順序が変わっているので、この場合には当てはまらない。しかしながら、分散・バラツキの場合は、順序が変わることは問題がなく、有用であることに変わりはない。

# 順序を保つ非負化関数

ドメイン(始域)が実数全体であって、コドメイン(終域)が非負であるような関数の例としては、
指数関数やアークタンジェント関数を用いたものがある。

以下では$g(x) = arctan(x) + \pi/2$という関数を例として用いる。

ドメインが実数全体である分布として正規分布をもう一度振り返る。

```{r}
a <- -3
b <- 2
x <- seq(from=-20,to=20,length=1000)
y <- dnorm(x,a,b)
plot(x,y,type="l")
```

閾値を描きいれ、理論的な$\mu$値を計算する。
```{r}
thres <- 0.8
plot(x,y,type="l")
abline(v=thres,col=2)
p0.8 <- pnorm(thres,a,b)
1-p0.8
```
累積分布を描く。

```{r}
t <- seq(from=-10,to=10,length=1000)
z <- pnorm(t,a,b)
plot(t,z,type="l",ylim=c(0,1))
abline(v=thres,col=2)
abline(h=p0.8,col=3)
```

この累積分布グラフの縦軸と横軸を入れ替えると次のようになる。

横軸はベータ分布のときと同じで0-1の範囲であるが、縦軸は、実数全体になる。

```{r}
plot(z,t,type="l")
abline(h=thres,col=2)
abline(v=p0.8,col=3)
```

階段状の関数が同じように描ける。

```{r}
w <- rep(0,length(z))
w[which(t>=thres)] = 1
matplot(z,cbind(t,w),type="l",col=c(1,4))
abline(h=thres,col=2)
abline(v=p0.8,col=3)
```

ここで黒い曲線が負の値を取っているのが不都合なのであった。

この縦軸の値を$g(u) = arctan(u)+\frac{\pi}{2}$で変換する。

閾値もgで変換する。
```{r}
g <- function(x){atan(x)+pi/2}
w <- rep(0,length(z))
w[which(t>=thres)] = 1
v <- g(t)
matplot(z,cbind(t,w,v),type="l",col=c(1,4,6))
abline(h=thres,col=2)
abline(v=p0.8,col=3)
abline(h=g(thres),col=2)
```

知りたい値に対応する長方形を描けば次のようになる。

```{r}
plot(z,v,type="l")
abline(v=p0.8,col=3)
abline(h=g(thres),col=2)
rect(p0.8,0,1,g(thres),col=6)
rect(p0.8,0,1,1,col=5)
```

上限を定めるのはここに描かれている黒い曲線の下面積を変換した閾値で割った値である。

Rの数値積分関数で下面積を求め、それを理論値と比較する。
積分したい範囲は、横軸がpnorm(t,a,b)、縦軸が、g(t)であるので、以下のコードでは、数値積分関数integrate()を使うにあたり、pnorm()の逆関数であるqnorm()を使って、『このグラフでの』横軸の値が、『このグラフでの』縦軸の値を返すような関数gf()を定義して、それを数値積分している。
```{r}
a <- -3
b <- 2
gf <- function(t){
  g(qnorm(t,a,b))
  
}
plot(z,gf(z))
area.below <- integrate(gf,0,1)[[1]]
area.below/gf(thres)
1-p0.8 <= area.below/gf(thres) 
```

閾値を本来の正規分布のドメイン全体でいろいろに変えれば、真の確率が、関数gを用いた上限によって抑えられる様子が描ける。

```{r}
up.bd <- area.below/g(t)
par(mfcol=c(1,2))
matplot(t,cbind(1-z,up.bd),type="l")
matplot(t,cbind(1-z,up.bd),type="l",ylim=c(-1.1,1.1))
par(mfcol=c(1,1))
```
ここで、以下の図を見直そう。
```{r}
plot(z,v,type="l")
abline(v=p0.8,col=3)
abline(h=g(thres),col=2)
rect(p0.8,0,1,g(thres),col=6)
rect(p0.8,0,1,1,col=5)
```
確かに紫の長方形の面積は黒線の下面積よりは小さいが、緑の縦線より右側のみの積分をしても、やはり紫の面積の上限を与えるし、下面積の全体よりも真値に近い。

それを示すのが以下の図(あまりに小さい閾値だと、数値積分の計算限界の兼ね合いで、不等式が満足されないが、ほどほどのところまでならば、不等式が満足されていることが示されている)。

以下のコードでは、積分計算に入れない範囲で黒線の縦軸の値が0となるような関数gf.()を定義してその数値積分を行っている。

```{r}
t <- seq(from=-1,to=2,length=100)
z <- pnorm(t,a,b)
up.bd <- area.below/g(t)
up.bd.2 <- rep(NA,length(t))
for(i in 1:length(up.bd.2)){
  k <- g(t[i])
  gf. <- function(x){
    tmp <- gf(x)
    tmp[which(tmp<k)] <- 0
    tmp
  }
  area.below.2 <- integrate(gf.,0,1,subdivisions=100)[[1]]
  up.bd.2[i] <- area.below.2/g(t[i])
}
matplot(t,cbind(1-z,up.bd,up.bd.2),type="l")
matplot(t,cbind(1-z,up.bd,up.bd.2),type="l",ylim=c(0,0.1))
matplot(t,cbind(1-z,up.bd,up.bd.2),type="l",ylim=c(0,0.01))
```

今、「使うとよいところ」だけを使って黒線の下面積を計算した。
このやり方の見方を変えると、「使わずに済ませたかったところ」の黒線が、グラフで描かれたカーブのところではなくて、値0にあるものとして積分した、ということである。

したがって、黒線に変換するときに用いたアークタンジェントを使った関数の代りに、一部はこの変換そのものを、他の部分は0となるような不連続な変換を用いても、Chebyshevの不等式は成り立つことがわかる。
ただし、その変換ではf(x)の大小順序を守るようなgであること、g(f(x))は非負であることは、処理の手続きから考えて必要である。

ちなみに、$|f(x)-E|$としたときのgは"平均からの距離”という関数だし、$(f(x)-E)^2$として分散を使えるようにしたときのgは$g(v)=v^2$という関数である。

また、Markovの不等式は、変換gをする必要がなく、正の範囲にある関数fについて、$g(v)=v$という場合に相当する。

それを意識して、Cehbyshevの不等式を表すと、

$$
\mu(\{x \in X ; f(x) \ge t\}) \le \frac{1}{g(t)}\int_{X} g(f(x)) dx
$$
となる。

ここで、

この式では、xに1次元数直線を考えて重み付き積分をしているが、もっと一般的にして、面積にしたり、体積にしたり、積分はできるがきれいな空間ではないような場合でも同じことなので、
$$
\mu(\{x \in X ; f(x) \ge t\}) \le \frac{1}{g(t)}\int_{X} g(f(x)) d\mu
$$
のように書くこともできる。

また、黒面積として使いたいところだけを使う、という条件を関数gの場合分けで書かずに、以下のように、積分範囲に制限を加えて表現することもできる。

$$
\mu(\{x \in X ; f(x) \ge t\}) \le \frac{1}{g(t)}\int_{f(x) \ge t} g(f(x)) d\mu
$$
または
$$
\mu(\{x \in X ; f(x) \ge t\}) \le \frac{1}{g(t)}\int_{g(f(x))\ge g(t)} g(f(x)) d\mu
$$