趣味の研究

趣味の数学を公開しています。初めての方は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