チェビシェフ不等式で遊ぶ
arXivに投稿したプレプリントの内容を確認するためのpythonサンプルコードです。
まずは、結果のグラフから。
このグラフは、平均値からの距離が以上の領域における確率密度関数の最大値が既知の場合、その領域全体の確率の上限を算出した結果です。
(この結果を1から引けば、]に含まれる確率の下限が求まります)
"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