趣味の研究

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

The graph of probability bound

I show the graph of probability bound.

Theorem1

Let X \in \mathcal{R} be a continuous random variable with expected value \mu and variance \sigma^2.

Let f(X) be a probability density function.

Let m_k be m_k=\sup_{|x-\mu|\geq k\sigma} f(x).

Then, 

Pr(|x-\mu|\geq k\sigma)\leq\alpha (m_k\sigma)

where

\alpha is a real root of 3-order equation.

N(x)=\frac{1}{2\pi}x^3+k x^2+e k^2 x -\frac{e}{m_k\sigma}=0.

 

Corollary 1.

Let X \in \mathcal{Z} be a discrete random variable with expected value \mu and variance \sigma^2.

Let p(X) be a probability mass function.

Let m_k be m_k=\sup_{\{x\leq\lceil {\mu-k\sigma} \rceil\}\cup\{x\geq\lfloor {\mu+k\sigma} \rfloor\}} p(x).

Then, 

Pr(|x-\mu|\geq k\sigma)\leq\alpha (m_kr\sigma)+p(\lceil\mu-k\sigma\rceil)

where

\alpha is a real root of 3-order equation.

N(x)=\frac{1}{2\pi}x^3+\frac{k}{r} x^2+e{(\frac{k}{r})}^2 x -\frac{e}{m_kr\sigma}=0.

r=\sqrt{1+\frac{1}{12\sigma^2}}.

 

About the proof, please see the article on April 11th.

The example of normal distribution.

f:id:hobbymath:20180414214851p:plain

The example of Poisson distribution.

f:id:hobbymath:20180414221908p:plain

"Probability" is the result of actual normal distribution, "Chebyshev" is Chebyshev inequality, and "New bound" is the result of Theorem 1. and Corollary 1.

 

The python code is 

 

# -*- coding: utf-8 -*-
"""
Created on Tue Apr 10 22:09:53 2018

@author: HobbyMath
"""

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

k_max =50
sigma=np.sqrt(3)
mu=3
Pr = np.ones(k_max , dtype=float)
chebyshev= np.ones(k_max , dtype=float)
bound= np.ones(k_max , dtype=float)

t = sympy.Symbol('t')
for i in range(1,k_max + 1):
print(i)
delta = 0.1
k = i * delta + 1
pos = k * sigma + mu
pos2 = - k* sigma + mu

pos_int = np.floor(k * sigma + mu) + 1
pos2_int = np.floor(-k * sigma + mu)

##########normal distribution######
# m = stats.norm.pdf(x=pos, loc=mu, scale=sigma)
# Pr[i] = 1 - stats.norm.cdf(x=pos, loc=mu, scale=sigma) + stats.norm.cdf(x=pos2, loc=mu, scale=sigma)

#########Poisson distribution######
m = max(stats.poisson.pmf(pos_int, mu), stats.poisson.pmf(pos2_int - 1, mu))
Pr[i] = 1 - stats.poisson.cdf(pos, mu) + stats.poisson.cdf(pos2, mu)

#########binomal distribution######
# n=10
# p=0.5
# sigma = np.sqrt(n * p * (1 - p))
# mu=n * p
# m=max(stats.binom.pmf(pos_int, n, p), stats.binom.pmf(pos2_int - 1, n, p))
# Pr[i] = 1 - stats.binom.cdf(pos, n, p) + stats.binom.cdf(pos2, n, p)

chebyshev[i] = 1 / k**2
r = np.sqrt(1 + 1 / (12 * sigma ** 2))
#r = 1
eq = 0.5 / np.pi * t**3 + k * t**2 / r + np.e * k**2 / r**2 * t - np.e / (m * r * sigma)
tmp = np.array(sympy.solve(eq), dtype=complex)
# bound[i]=np.abs(tmp[0]) * m * r * sigma
bound[i]=np.abs(tmp[0]) * m * r * sigma + stats.poisson.pmf(pos2_int , mu)


plt.title("Poisson distribution(mu=1, sigma=1")
plt.ylabel("probability")
plt.xlabel("k")
#plt.yscale('log')

prange = np.arange(1 ,k_max)
x = prange.astype(np.float64) * delta + 1.0

p1 = plt.plot(x, Pr[prange])
p2 = plt.plot(x, chebyshev[prange])

p3 = plt.plot(x, bound[prange])

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