import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
from math import atan

def gaussian_pdf(x,mu,sigma):
    return(1/(np.sqrt(2*np.pi)*sigma)*np.exp(-(x-mu)**2/(2*sigma**2)))

def lorentzian_pdf(x,x0,w):
    return(1/np.pi*w/2/((w/2)**2+(x-x0)**2))

def gaussian_cdf(x,mu,sigma):
    return(0.5*(1+erf((x-mu)/(sigma*np.sqrt(2)))))

def lorentzian_cdf(x,x0,w):
    return(1/np.pi*(np.arctan(2*(x-x0)/w)+np.pi/2))

mu=0
x0=0
w=1
sigma=1/(2*np.sqrt(2*np.log(2)))

x=np.arange(mu-5*w,mu+5*w,0.01)
plt.plot(x,gaussian_pdf(x,mu,sigma),c='red',label='Gaussian')
plt.plot(x,lorentzian_pdf(x,x0,w),label='Lorentzian')
plt.xlim(mu-5*w,mu+5*w)
plt.xlabel('x')
plt.ylabel('pdf(x)')
plt.legend()
plt.show()

plt.plot(x,gaussian_cdf(x,mu,sigma),c='red', label='Gaussian')
plt.plot(x,lorentzian_cdf(x,x0,w),label='Lorentzian')
plt.xlim(mu-5*w,mu+5*w)
plt.ylim(0,)
plt.xlabel('x')
plt.ylabel('cdf(x)')
plt.legend()
plt.show()

print('Gaussian: P(|x|>2)   = ',1-gaussian_cdf(2,mu,sigma))
print('Lorentzian: P(|x|>2) = ',1-lorentzian_cdf(2,x0,w))      

Gaussian: P(|x|>2)   =  1.240773008670537e-06
Lorentzian: P(|x|>2) =  0.07797913037736937