In [1]:
# w00-in class code
import warnings
import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.stats import expon
import random
import os
# Normal distribution
# P(x) = 1/sqrt(2*pi*sigma^2) * exp[ -(x-mu)^2/(2*sigma^2) ]
#
mu = 5
sigma = 2
Nsamples = 100000
# rsv=random variable sample
# xs = samples of normal dist
xs = scipy.stats.norm.rvs(mu, sigma, Nsamples)
# calculate the cdf of the xs samples
ys = scipy.stats.norm.cdf(xs, mu, sigma)
fig, axes = plt.subplots(1,2, figsize=(9, 3))
axes[0].hist(xs, bins=50)
axes[0].set_title("X Samples")
axes[1].hist(ys, bins=50)
axes[1].set_title("Fx Samples")
Out[1]:
Text(0.5, 1.0, 'Fx Samples')
In [2]:
ld = 0.2
# rsv=random variable sample
# xs = samples of normal dist
xs = scipy.stats.expon.rvs(scale = 1/ld, size=Nsamples)
# calculate the cdf of the xs samples
ys = scipy.stats.expon.cdf(xs, scale = 1/ld)
fig, axes = plt.subplots(1,2, figsize=(9, 3))
axes[0].hist(xs, bins=50)
axes[0].set_title("X Samples")
axes[1].hist(ys, bins=50)
axes[1].set_title("Fx Samples")
Out[2]:
Text(0.5, 1.0, 'Fx Samples')
In [3]:
# w00-in class code
import warnings
import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.stats import expon
import random
import os
Nsamples = 4000
fig, axes = plt.subplots(1,3, figsize=(9, 3))
# Normal distribution
# P(x) = 1/sqrt(2*pi*sigma^2) * exp[ -(x-mu)^2/(2*sigma^2) ]
#
mu = 5
sigma = 2
xs = scipy.stats.norm.rvs(mu, sigma, Nsamples)
ys = scipy.stats.norm.cdf(xs, mu, sigma)
axes[0].hist(xs, bins=50)
axes[0].set_title("X Samples")
axes[1].hist(xs, density=True, cumulative=True, bins=50)
axes[1].hist(xs, density=True, bins=50)
axes[1].set_title("X PDF/CDF")
x = np.linspace(norm.ppf(0.001, loc=mu, scale=sigma),norm.ppf(0.99, loc=mu, scale=sigma), 100)
axes[1].plot(x, norm.pdf(x, loc=mu, scale=sigma),'g-', lw=3, alpha=0.6)
axes[1].plot(x, norm.cdf(x, loc=mu, scale=sigma),'g-', lw=3, alpha=0.6)
axes[2].hist(ys, bins=50)
axes[2].set_title("CDF=F_X(X) samples")
plt.show()
plt.savefig('Normal_CDFisUniform.png')
plt.close(fig)
<Figure size 432x288 with 0 Axes>
In [ ]:
In [4]:
fig, axes = plt.subplots(1,3, figsize=(9, 3))
Nsamples = 10000000
# Exponential distribution
# P(x) = ld * exp(-ld*x)
#
ld = 1.3
xs = scipy.stats.expon.rvs(scale = 1/ld, size=Nsamples)
ys = scipy.stats.expon.cdf(xs, scale = 1/ld)
axes[0].hist(xs, bins=50)
axes[0].set_title("X Samples")
axes[1].hist(xs, density=True, cumulative=True, bins=50)
axes[1].hist(xs, density=True, bins=50)
axes[1].set_title("X PDF/CDF")
x = np.linspace(expon.ppf(0.001, scale=1/ld),expon.ppf(0.99, scale=1/ld), 100)
axes[1].plot(x, expon.pdf(x, scale=1/ld),'g-', lw=3, alpha=0.6)
axes[1].plot(x, expon.cdf(x, scale=1/ld),'g-', lw=3, alpha=0.6)
axes[2].hist(ys, bins=50)
axes[2].set_title("CDF=F_X(X) samples")
plt.show()
plt.savefig('Exponential_CDFisUniform.png')
plt.close(fig)
<Figure size 432x288 with 0 Axes>
In [5]:
# Calculate CDF from PDF numerically
#
fig, axes = plt.subplots(2,3, figsize=(20, 10))
# inverse sampling for Normal distribution N(0,1)
x = np.linspace(norm.ppf(0.001),norm.ppf(0.999), 100)
pdf = norm.pdf(x)
print(x)
print(pdf)
Nbin = len(pdf)
bsize = (x[Nbin-1] - x[0])/Nbin
print("xmin", x[0], "xmax", x[Nbin-1], "Nbin", Nbin, "bsize", bsize)
axes[0][0].plot(pdf)
axes[0][0].set_title("Nomal PDF")
[-3.09023231 -3.02780337 -2.96537444 -2.9029455 -2.84051656 -2.77808763 -2.71565869 -2.65322976 -2.59080082 -2.52837189 -2.46594295 -2.40351402 -2.34108508 -2.27865614 -2.21622721 -2.15379827 -2.09136934 -2.0289404 -1.96651147 -1.90408253 -1.8416536 -1.77922466 -1.71679573 -1.65436679 -1.59193785 -1.52950892 -1.46707998 -1.40465105 -1.34222211 -1.27979318 -1.21736424 -1.15493531 -1.09250637 -1.03007744 -0.9676485 -0.90521956 -0.84279063 -0.78036169 -0.71793276 -0.65550382 -0.59307489 -0.53064595 -0.46821702 -0.40578808 -0.34335915 -0.28093021 -0.21850127 -0.15607234 -0.0936434 -0.03121447 0.03121447 0.0936434 0.15607234 0.21850127 0.28093021 0.34335915 0.40578808 0.46821702 0.53064595 0.59307489 0.65550382 0.71793276 0.78036169 0.84279063 0.90521956 0.9676485 1.03007744 1.09250637 1.15493531 1.21736424 1.27979318 1.34222211 1.40465105 1.46707998 1.52950892 1.59193785 1.65436679 1.71679573 1.77922466 1.8416536 1.90408253 1.96651147 2.0289404 2.09136934 2.15379827 2.21622721 2.27865614 2.34108508 2.40351402 2.46594295 2.52837189 2.59080082 2.65322976 2.71565869 2.77808763 2.84051656 2.9029455 2.96537444 3.02780337 3.09023231] [0.00336709 0.00407561 0.00491403 0.00590188 0.00706074 0.00841429 0.00998831 0.01181066 0.01391117 0.01632152 0.01907502 0.02220632 0.02575109 0.02974556 0.03422599 0.0392281 0.04478637 0.05093331 0.05769861 0.06510827 0.0731837 0.08194076 0.0913888 0.10152977 0.11235727 0.1238558 0.136 0.14875408 0.16207135 0.17589399 0.19015298 0.20476827 0.21964919 0.23469504 0.24979609 0.26483461 0.27968633 0.29422199 0.30830915 0.32181412 0.33460403 0.34654899 0.35752424 0.36741234 0.37610524 0.38350623 0.38953174 0.39411292 0.39719693 0.39874797 0.39874797 0.39719693 0.39411292 0.38953174 0.38350623 0.37610524 0.36741234 0.35752424 0.34654899 0.33460403 0.32181412 0.30830915 0.29422199 0.27968633 0.26483461 0.24979609 0.23469504 0.21964919 0.20476827 0.19015298 0.17589399 0.16207135 0.14875408 0.136 0.1238558 0.11235727 0.10152977 0.0913888 0.08194076 0.0731837 0.06510827 0.05769861 0.05093331 0.04478637 0.0392281 0.03422599 0.02974556 0.02575109 0.02220632 0.01907502 0.01632152 0.01391117 0.01181066 0.00998831 0.00841429 0.00706074 0.00590188 0.00491403 0.00407561 0.00336709] xmin -3.090232306167813 xmax 3.090232306167813 Nbin 100 bsize 0.06180464612335626
Out[5]:
Text(0.5, 1.0, 'Nomal PDF')
In [6]:
fig, axes = plt.subplots(2,3, figsize=(20, 10))
# inverse sampling
#
# We need the CDF
x = np.linspace(norm.ppf(0.001),norm.ppf(0.999), 100)
pdf = norm.pdf(x)
#print(pdf)
Nbin = len(pdf)
bsize = (x[Nbin-1] - x[0])/Nbin
print("xmin", x[0], "xmax", x[Nbin-1], "Nbin", Nbin, "bsize", bsize)
axes[0][0].plot(pdf)
axes[0][0].set_title("Nomal PDF")
# CDF (not analytically)
#
# CDF(x) = CDF(x-b) + PDF(x) * b
#
cdf = [None]*Nbin
cdf[0] = 0
for b in range(1,Nbin-1):
cdf[b] = cdf[b-1] + pdf[b]*bsize
axes[0][1].plot(cdf)
axes[0][1].set_title("Nomal CDF")
axes[0][2].plot(pdf)
axes[0][2].set_title("Nomal PDF")
axes[0][2].plot(cdf)
axes[0][2].set_title("Nomal CDF")
xmin -3.090232306167813 xmax 3.090232306167813 Nbin 100 bsize 0.06180464612335626
Out[6]:
Text(0.5, 1.0, 'Nomal CDF')
In [ ]:
In [ ]:
In [ ]: