趣味の研究

趣味の数学を公開しています。初めての方はaboutをご覧ください。

チェビシェフ不等式で遊ぶ

arXivに投稿したプレプリントの内容を確認するためのpythonサンプルコードです。

まずは、結果のグラフから。

f:id:hobbymath:20180907223409p:plain

このグラフは、平均値\muからの距離が\epsilon\sigma以上の領域における確率密度関数の最大値が既知の場合、その領域全体の確率の上限を算出した結果です。

(この結果を1から引けば、[\mu-\epsilon\sigma, \mu+\epsilon\sigma]に含まれる確率の下限が求まります)

"Probability"は正規分布の場合に計算した例、"Chebyshev"がチェビシェフ不等式に基づいて算出した確率の上限、"New bound"が今回のプレプリントの内容をもとに算出した確率の上限です。

 

 ここからpythonのサンプルコードです。

 

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 30 10:04:20 2018

@author: Nishiyama
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import sympy

k_max =50
delta = 0.1

sigma = 1 #standard deviation
mu= 1 #mean
Probability = np.zeros(k_max , dtype=float)
Chebyshev = np.zeros(k_max , dtype=float)
bound= np.zeros(k_max , dtype=float)
epsilon = np.zeros(k_max , dtype=float)


########Calculate probability bound##########
def calc_probability(mu, sigma, epsilon, sup):
#sup : #supremum of PDF

t = sympy.Symbol('t')
m_e = sigma * sup
eq = 0.5 / (np.pi * np.e) * t**3 + epsilon / np.e * t**2 + epsilon**2 * t - 1/ m_e
tmp = np.array(sympy.solve(eq), dtype=complex)
return np.abs(tmp[0]) * m_e

#######Example of normal distribution#######
for i in range(0,k_max):
print(i)

epsilon[i] = i * delta + 1
pos_R = epsilon[i] * sigma + mu
pos_L = -epsilon[i]* sigma + mu

#calculate supremum of PDF
sup= max(stats.norm.pdf(x=pos_L, loc=mu, scale=sigma), stats.norm.pdf(x=pos_R, loc=mu, scale=sigma))
tmp_probability = stats.norm.cdf(x=pos_R, loc=mu, scale=sigma) - stats.norm.cdf(x=pos_L, loc=mu, scale=sigma)
Probability[i] = 1 - tmp_probability
bound[i]= calc_probability(mu, sigma, epsilon[i], sup)

Chebyshev = 1 / epsilon**2

plt.show()
plt.ylabel("probability")
plt.xlabel("epsilon")

prange = np.arange(0 ,k_max)

p1 = plt.plot(epsilon[prange], Probability[prange])
p2 = plt.plot(epsilon[prange], Chebyshev[prange])
p3 = plt.plot(epsilon[prange], bound[prange])

plt.legend((p1[0], p2[0], p3[0]), ("Probability", "Chebyshev", "New bound"), loc=5)

 

プレプリントのリンク(月曜日にグラフを微修正したバージョンをアップロード予定です)

[1808.10770] Improved Chebyshev inequality: new probability bounds with known supremum of PDF

 

チェビシェフ不等式に関する結果のまとめ

チェビシェフ不等式に関する考察のまとめをプレプリントとして投稿しました。


概要

平均値より一定以上離れた集合における確率密度関数の上限(最大値)が既知の場合、確率に対してより厳しい不等式を導出することができる、というものです。

https://osf.io/h9zfn/

内容のまとめ

これまで書いた内容を順番にまとめておきます。

内容の整理は、ぼちぼちと進めます。

 

1. 確率微分方程式(ブラウン運動モデル)を用いたコラッツ問題の解析

  深堀したいのですが、コラッツ問題にこのモデルを適用する正当性がうまく

 説明できずにいます。

 

2.素数定理と情報エントロピー

 お遊び的な内容です。

 

3.Renyiエントロピー不等式の関数不等式への活用

 

4.(Renyi)エントロピーのentropy power inequalityに関する考察

 

5.Chebyshev不等式の拡張

 地味な内容ですが、こちらもそのうち深堀したいです。

 

6. wavelet変換と混合確率分布の関係に関する考察

 あまり厳密ではありません。

 

7. ダイバージェンスに関する考察(arXivに投稿しました)

[1808.06482] Divergence functions in dually flat spaces and their properties

[1808.06148] Generalized Bregman and Jensen divergences which include some f-divergences

 

 

カルバック・ライブラー情報量と平行四辺形

機械学習の分野でよく用いられるクロスエントロピーや、カルバック・ライブラー情報量に関する関係を導き出したので、要約を書いておきます。

今回見つけたことは、ざっくり言うと、

4つの確率分布が、「平行四辺形」の各頂点にあるとき、確率分布間の「距離」に対し、平行四辺形の法則(中線定理)が成り立つ

というものです。「」をつけて書いたのは、通常の意味の平行四辺形や、距離とは異なるからです。

なるべく式を使わずに内容を書いていきます。

詳細は6/21の記事(清書版)、6/22の記事をご覧ください。

 

確率分布 

\begin{align}p(x)=\exp(\sum_i \theta^iF_i(x)-\psi(\theta) )\tag{1}\end{align}

を考えます。ここで、\theta^iのiはベキ乗の意味ではなく、変数の添え字です。上についていることには意味がありますが、ここでは無視してください。

このような形の確率分布は、指数型分布族と呼ばれ、正規分布、指数分布、ポアソン分布など

多くの分布が該当します。

 

次にカルバック・ライブラー情報量を考えます。

これは、確率分布同士がどれくらい似ているかを表す量で、距離のようなものです。

これをd_{KL}(p||q)と書きます。

クロスエントロピーとの関係は、クロスエントロピーH(p,q), エントロピーH(p)と書くと、

 

d_{KL}(p||q)=H(p,q)-H(p)です。

この量は、p,qに対して対称ではないので、pとqをひっくり返したものを足して、あらたに

\begin{align}D_G(p,q)\equiv d_{KL}(p,q)+d_{KL}(q,p)\tag{2}\end{align}

という量を考えます。

 

次に、確率分布を4つ考え、p, q, r, s

とします。これらの確率分布は、すべて(1)式の形で書けるとし、パラメータ\theta\psi(\theta)の部分だけが異なるとします(つまりF_i(x)の部分は同じ)。

 

いよいよ今回見つけた関係式の説明に入ります。

上の確率分布の間に

\theta^i(p)+\theta^i(s)=\theta^i(q)+\theta^i(r)

または、

E[F_i]_p+E[F_i]_s=E[F_i]_q+E[F_i]_rの関係式がすべてのiに対して成り立つとします。F_i(x)xは省略して記載しました。

 

E[ ]_pは確率分布pにおける期待値の意味です。この関係式は、p, q, r, sが普通のベクトルの場合、p, q, r, sが平行四辺形の頂点であることを意味します。

 

もし、確率分布間に上のいずれかの関係が成り立つと、

(2)で定義したD_Gに対して

\begin{align}D_G(p,q)+D_G(q,s)+D_G(s,r)+D_G(r,p)=D_G(p,s)+D_G(q,r)\tag{3}\end{align}

が成り立ちます。

(3)、D_G(p,q)を距離の2乗とみなしたとき、平行四辺形p, q, r, sの辺の2乗の和が対角線の2乗の和と等しいことを意味しています。

これは、幾何で習う平行四辺形の法則(中線定理)の関係式そのものです。

確率分布でも同じ関係が成り立つのは、なんだか不思議です。

 

 

情報幾何学におけるダイバージェンス

この記事は、6/21に書いた記事を清書したものです。

 

◇概要

情報幾何学は、確率分布を統計多様体上の点とみなし、統計多様体の幾何的な性質を解析する学問である。統計多様体が双対平坦な場合には、確率分布間の「距離の2乗」に該当する正準ダイバージェンスを導入することができる。3つの確率分布および、正準ダイバージェンスに対し、余弦定理に対応する三角関係式が成り立つことが知られている。特殊な場合において、三角関係式は拡張ピタゴラスの定理に帰着し、射影に関して基礎的な定理となる。

本記事では、三角関係式を用いて、主に二つの性質を示す。

まず、正準ダイバージェンスから新しいダイバージェンス(幾何ダイバージェンス(仮称)を導入し、正準ダイバージェンス、幾何ダイバージェンスともに、2つのアファイン座標の測地線上で、単調性を有することを示す。

次に、確率分布同士の「ベクトル和」演算を導入し、4つの確率分布が「平行四辺形」の頂点上にあるとき、幾何ダイバージェンスに対して、平行四辺形の法則(中線定理)が成り立つことを示す。

◇基本的な定義[1]

\thetaを統計多様体におけるアファイン座標、\etaをその双対座標(期待値座標)とする。

統計多様体の二点p, qを考え、q\theta座標を\theta(p), q\eta座標を\eta(q)とする。

このとき、二点間の正準ダイバージェンス

D(p||q)\equiv \psi(\theta(p) )+\phi(\eta(q) )-\theta(p)^i\eta_i(q)

で定義する。

上付き添え字、下付き添え字がペアで現れたときは、\sumの記号を省略するものとする。(Einsteinの縮約規則)

\phi(\theta), \psi(\eta)はポテンシャルであり、互いにLegendre変換の関係にある。即ち、

\phi(\eta)=\theta^i\eta_i-\psi(\theta)

以下では、\psi(\theta(p) )\psi(p), \phi(\eta(p) )\phi(p)のように略記する。

\theta,\eta,\psi,\phiは次の関係式を満たす。

\theta^i=\partial^i\phi\eta_i=\partial_i\psig_{ij}=\partial_i\eta_jg^{ij}=\partial^i\theta^j

ここで、g_{ij}はFisher情報量から導出されるRiemann計量とする。

指数型分布族の場合、おおむね\psiは負の自由エネルギーを表し、\phiは負のエントロピーを表す。

詳細は参考文献[1]を参照のこと。

 

 ◇幾何ダイバージェンスの定義

D_{G}(p,q)\equiv D(p||q)+D(q||p)と定義し、D_G(p,q)を幾何ダイバージェンスと呼ぶことにする。

幾何ダイバージェンスは、Euclid空間のL2ノルムと類似した性質を示す(後述)。

定義から、幾何ダイバージェンスは点p,qに対し対称である。

また、指数型分布族に対しては、幾何ダイバージェンスは、Kullbalk-Leiblerダイバージェンスの対称和にほぼ一致する。

D_G(p,q)=d_{KL}(p||q)+d_{KL}(q||p)

 

命題1. 

幾何ダイバージェンスは、アファイン座標を用いてD_G(p,q)=(\eta_i(q)-\eta_i(p) )(\theta^i(q)-\theta^i(p) )と表される。

また、D_G(p,q)\geq 0であり、D_G(p,q)=0となるのは、p=qのみである。

 

証明 

正準ダイバージェンスの三角関係式D(p||q)+D(q||r)-D(p||r)=(\eta_i(r)-\eta_i(q) )(\theta^i(p)-\theta^i(q) )

において、p=rと置くことで、D_G(p,q)=(\eta_i(q)-\eta_i(p) )(\theta^i(q)-\theta^i(p) )を得る。

後半は、正準ダイバージェンスD(p||q)\geq 0(等号成立はp=qのみ)から従う。

 

Riemann空間における微小距離の2乗はd\eta_id\theta^iと書けることから、幾何ダイバージェンスも距離の2乗のような振る舞いをすると予想される。

 

◇測地線における単調性

正準ダイバージェンス、幾何ダイバージェンスとも測地線の方向ベクトルに対し、単調増加関数となる。これは、距離において望まれる性質である。

 

定理 1.

正準ダイバージェンスD(p||q)は、pを固定し、qを\etaもしくは、\theta-測地線に沿って動かした場合、測地線の方向ベクトルに対し、単調増加関数となる。

qを固定してpを動かした場合も同様である。

 

証明

\begin{align}\partial_{q,i}D(p||q)=g_{ij}(q)(\theta^j(q)-\theta^j(p) )\tag{1}\end{align}  

\begin{align}\partial_q^iD(p||q)=\theta^i(q)-\theta^i(p)\tag{2} \end{align}

を用いる。

\theta-測地線の場合、\theta^i(q)=\theta^i(p)+a^itと書けるので、

(1)を用いると、

\begin{align}\frac{dD(p||q)}{dt}=a^j\partial_{q,j}D(p||q)=ta^jg_{ij}(q)a^i\tag{3}\end{align} 

右辺は、計量の正定値性より、tの符号と一致する。よって、単調性が示せた。

\eta-測地線は、\eta_i(q)=\eta_i(p)+a_itと書ける。(2)より、

\begin{align}\frac{dD(p||q)}{dt}=a_i\partial_q^iD(p||q)=a_i(\theta^i(q)-\theta^i(p) )=\frac{1}{t}(\eta_i(q)-\eta_i(p)(\theta^i(q)-\theta^i(p) )=\frac{1}{t}D_G(p,q)\end{align}

\begin{align}\tag{4}\end{align}となる。

命題1より、D_G\geq 0であるから、右辺はtの符号と一致し、単調性が示せた。

 
pに関して
\begin{align}\partial_p^iD(p||q)=g^{ij}(p)(\eta_j(q)-\eta_j(p) )\tag{5}\end{align}  
\begin{align}\partial_{p,i}=\theta^i(q)-\theta^i(p)\tag{6}\end{align}      
となることから、同様に証明できる。
 
式(3)を\frac{dD(p||q)}{dt}=\frac{1}{t}\{g_{ij}(a^it)(a^jt)\}
と書き直すと、g_{ij}(a^it)(a^jt)はRiemann空間における距離の2乗のようなものなので、式(4)と対を成す。
 
系1.
幾何ダイバージェンスD_G(p,q)は、\theta, \eta-測地線に沿ってp,qを動かしたとき単調性を示す。
定理1より、即座に導かれる。
 

内積、ベクトル和の定義

内積を以下の式で定義する。
\langle q,r\rangle_p\equiv\frac{1}{2}(\eta_i(q)-\eta_i(p) )(\theta^i(r)-\theta^i(p) )+\{q\leftrightarrow r\}
\leftrightarrowは記号の入れ替えを意味する。
特に、q=rであれば、\langle q,q\rangle_p=D_G(p,q)である。
 
点q,rのベクトル和を以下の式で定義する。
s=(q+r)_\thetaと書いたとき、すべてのiに対し、\theta^i(s)=\theta^i(q)+\theta^i(r)を満たす点であるとする。
同様に、s=(q+r)_\etaと書いたとき、\eta_i(s)=\eta_i(q)+\eta_i(r)を満たす点であるとする。
このように導入された内積は、内積の線形の定義を満たさないことに注意。
系2. (余弦定理の拡張)
幾何ダイバージェンスは、Euclid空間のL2ノルムと同様、
\begin{align}D_G(p,q)+D_G(p,r)-2\langle q,r\rangle_p=D_G(q,r)\tag{7}\end{align}      
を満たす。
 
証明
正準ダイバージェンスの三角関係式そのものである。
 
定理 2.(ベクトル和の公式)
点p,q,r,sが(s+p)_\theta=(q+r)_\thetaもしくは、(s+p)_\eta=(q+r)_\etaを満たすとき、
\begin{align}D_G(q,s)+D_G(r,s)+2\langle q,r\rangle_p=D_G(p,s)\tag{8}\end{align}が成り立つ。  
 また、対称性から、pとsを入れ替えた式も成り立つ。
 
証明
(s+p)_\theta=(q+r)_\thetaの場合について証明する。
 
D_G(p,s)=(\theta^i(s)-\theta^i(p) )(\eta_i(s)-\eta_i(p) )=(\theta^i(r)-\theta^i(p) )(\eta_i(s)-\eta_i(p) ) + (\theta^i(q)-\theta^i(p) )(\eta_i(s)-\eta_i(p) )
D_G(q,s)=(\theta^i(s)-\theta^i(q) )(\eta_i(s)-\eta_i(q) )=(\theta^i(r)-\theta^i(p) )(\eta_i(s)-\eta_i(q) )
D_G(r,s)=(\theta^i(s)-\theta^i(r) )(\eta_i(s)-\eta_i(r) )=(\theta^i(q)-\theta^i(p) )(\eta_i(s)-\eta_i(r) )
 
を代入することにより、示すことができる。
この定理はp,q,r,sが平行四辺形の頂点である場合に対応する。
 
系3.(平行四辺形の定理の拡張)
点p,q,r,sが
(p+s)_\theta=(q+r)_\thetaもしくは、(p+s)_\eta=(q+r)_\etaを満たすとき、
\begin{align}D_G(p,q)+D_G(q,s)+D_G(s,r)+D_G(r,p)=D_G(q,r)+D_G(p,s)\tag{9}\end{align}が成り立つ。
 
これは、Euclid空間において、四角形p,q,r,sが平行四辺形の時に成り立つ平行四辺形の定理の拡張である。
 
証明
(8)に(7)を代入すればよい。
 
特に、p=sとおくと、
\theta(p)=\frac{1}{2}(\theta(q)+\theta(r) )もしくは、\eta(p)=\frac{1}{2}(\eta(q)+\eta(r) )を満たす点p,q,rに対し
D_G(q,p)+D_G(p,r)=\frac{1}{2}D_G(q,r)が成り立つ。
 

References

[1]甘利 俊一, "情報幾何とその応用"

Mixture distribution(e.g. GMM) and wavelet transform

Can we represent the arbitrary function by mixture distribution?

I thought this problem by using the analogy of wavelet transform.


Definition.

Let p(x) be a probability density  function in  \mathrm{R^n}.

p_{a,b}(x)=\Pi_k|a_k|^{-1}p(\frac{x-b}{a})

for distribution function p(x) 

p(x) indicates p(x_1,x_2,\cdot\cdot\cdot,x_n).

In the following, we omitte the arguments of function in the same way.

We confirm p_{a,b} is probability density  function easily.


Proposition

Let p(x)\in L^2(\mathrm{R^n}) be probability density  function.

and C^2 class function such that

C_p=\int\Pi_k|\omega_k| |\hat{p}(\omega)|^2d^n\omega<+\infty


We can expand any probability density  function f(x)\in L^2(\mathrm{R^n}),

f(x)=\int_{\mathrm{R^n}}d^na\int_{\mathrm{R^n}}d^nbF(a,b)p_{a,b}(x)     (1)


We represent the expansion coefficients as below.

F(a,b)=-\frac{\Pi_k|a_k|}{C_p}\int[\Pi_k\frac{d^2}{dt_k^2}f(t)] p_{a,b}(t)d^nt        (2)

We ommitte the integral range in the case the integral range is \mathrm{R^n}.

Where, \hat{p}(\omega) is Fourier transform of p(x).

We have Gaussian mixture model if p(x) is normal distribution p.d.f.

 

Proof.

We can prove in the same way as continuous wavelet transform.([1]Theorem4.4)

The right integral 

q(x)=\int d^na\int d^nbF(a,b)p_{a,b}(x) in (1) can be rewritten as the sum of convolution.

q(x)=\int d^na\int d^nbF(a,b)\Pi_k |a_k|^{-1} p(\frac{x-b}{a})=\int d^na\Pi_k |a_k|^{-1} F(a,\cdot\cdot)\star p_a(x)     

\star indicates convolution, the "." indicates the variable over which the convolution is caliculated.

We define p_a(x)=p(\frac{x}{a}).

In the same way, we have

F(a,b)=-\frac{1}{C_p}[\Pi_k\frac{d^2}{dt_k^2}f]\star \tilde{p_a}(b)

 where, \tilde{p_a}(x)=p_a(-x).
Combining the results, 
q(x)=-\frac{1}{C_p}\int d^na \Pi_k|a_k|^{-1}[\Pi_k\frac{d^2}{dt_k^2}f]\star\tilde{p_a}\star p_a(x).
Using Fourier transform formula, we have
\hat{p_a}(\omega)=\Pi_k a_k\hat{p}(a\omega)
and 
\hat{\tilde{p_a}}(\omega)=\hat{p_a}(-\omega)=\overline{\hat{p_a}(\omega)}.
\overline{p} indicates the complex conjugate of p and we use p_a is a real function.
Fourier transform of q(x) is
\hat{q}(\omega)=\frac{\hat{f}(\omega)}{C_p}\int d^na \Pi_k|a_k|^{-1}a_k^2\omega_k^2 |\hat{p}(a\omega)|^2.

The change of variable a'_k=\omega_k a_k and using |\hat{p}(\omega)|^2=|\hat{p}(-\omega)|^2, we have
\hat{q}(\omega)=\hat{f}(\omega).        (3)
By inverse Fourier transform to (3), we prove the Proposition.

Reference.
[1] St´ephane Mallat., "A Wavelet Tour of Signal Processing"

 

 

 

混合確率分布とウェーブレット変換

任意の確率密度関数は別の確率密度関数の混合確率分布で表されるのか?


という問題に対し、ウェーブレット変換の手法を用いて考えてみます。

 

⚪︎定義

\mathrm{R^n}で定義された確率密度関数p(x)に対し、

p_{a,b}(x)=\Pi_k|a_k|^{-1}p(\frac{x-b}{a})

と定義する。

また、p(x)は、p(x_1,x_2,\cdot\cdot\cdot,x_n)を、(\frac{x-b}{a})は、(\frac{x_1-b_1}{a_1},\frac{x_2-b_2}{a_2},\cdot\cdot\cdot,\frac{x_n-b_n}{a_n})を表すものとする。

p_{a,b}確率密度関数になることは容易に確かめることができる。

以下では、細字の添え字についてもベクトルを表すものとする。

 

○命題?

p(x),f(x)\mathrm{R^n}上で定義された二乗可積分確率密度関数とする。さらに、f(x)C^2級であるとする。

このとき、f(x)

f(x)=\int_{\mathrm{R^n}}d^na\int_{\mathrm{R^n}}d^nbF(a,b)p_{a,b}(x)

と展開可能である。

以下では、特に断りがない限り、積分範囲は\mathrm{R^n}であるとする。

展開係数は例えば、

F(a,b)=-\frac{\Pi_k|a_k|}{C_p}\int[\Pi_k\frac{d^2}{dt_k^2}f(t)] p_{a,b}(t)d^nt

C_p=\int\Pi_k|\omega_k| |\hat{p}(\omega)|^2d^n\omega

と表すことができる。

ここで、\hat{p}(\omega)p(x)フーリエ変換である。

特に、p(x)として、正規分布をとれば、混合ガウス分布の展開公式になる。

 

証明

連続ウェーブレット変換を導く場合([1]のTheorem4.4)と同様の考え方である。

q(x)=\int d^na\int d^nbF(a,b)p_{a,b}(x)

とおくと、

q(x)=\int d^na\int d^nbF(a,b)\Pi_k |a_k|^{-1} p(\frac{x-b}{a})=\int d^na\Pi_k |a_k|^{-1} F(a,\cdot\cdot)\star p_a(x)     

ここで、\starは畳み込みを表し、F(a,\cdot\cdot)内の\cdotは畳み込みに関する変数を表す。また、p_a(x)=p(\frac{x}{a})である。

同様に、

F(a,b)=-\frac{1}{C_p}[\Pi_k\frac{d^2}{dt_k^2}f]\star \tilde{p_a}(b)

 ここで、\tilde{p_a}(x)=p_a(-x)である。
合わせて、
q(x)=-\frac{1}{C_p}\int d^na \Pi_k|a_k|^{-1}[\Pi_k\frac{d^2}{dt_k^2}f]\star\tilde{p_a}\star p_a(x)
\hat{p_a}(\omega)=\Pi_k a_k\hat{p}(a\omega)
および、
\hat{\tilde{p_a}}(\omega)=\hat{p_a}(-\omega)=\overline{\hat{p_a}(\omega)}
を用いる。\overline{p}は、p複素共役である。最後の等式は、p_aが実数値関数であることから従う。
q(x)フーリエ変換すると、
\hat{q}(\omega)=\frac{\hat{f}(\omega)}{C_p}\int d^na \Pi_k|a_k|^{-1}a_k^2\omega_k^2 |\hat{p}(a\omega)|^2
a'_k=\omega_k |a_k|と変数変換して、
\hat{q}(\omega)=\hat{f}(\omega)
となるので、フーリエ逆変換して結果を得る。

追加:
Fourier変換の際、\omegaa積分を交換してますが、厳密には、まずa_kに関する積分区間を有限にします。
Youngの不等式により、\omegaに関する積分結果が有限になることを示せるので、Fubiniの定理により、積分順序を交換します。
そのあと、a_k積分範囲を\mathrm{R^n}にする極限をとることで、積分順序が交換可能であることを示します。

Reference.
[1] St´ephane Mallat., "A Wavelet Tour of Signal Processing"