Skip to content

Instantly share code, notes, and snippets.

@ofey404
Created October 16, 2019 02:52
Show Gist options
  • Save ofey404/b9513b152820534069217678a2c9fd4e to your computer and use it in GitHub Desktop.
Save ofey404/b9513b152820534069217678a2c9fd4e to your computer and use it in GitHub Desktop.
import numpy as np
from numpy.random import poisson
import matplotlib.pyplot as plt
import scipy.misc
import scipy.special
# Poisson distribution
# a discrete probability distribution that expresses the probability of a given
# number of events occurring in a fixed interval of time or space if these
# events occur with a known constant rate and independently of the time since
# the last event.
# Formula:
# P(k event in the interval) = (lambda^k * exp(-lambda)) / (k!)
# lambda: The average number of events in an interval
ax1 = plt.subplot(2, 2, 1)
lbd = 5
s = np.random.poisson(lbd, 100000)
count, bins, ignored = ax1.hist(s, 14, normed=True)
ax1.plot(bins, lbd**bins*np.exp(-lbd)/scipy.misc.factorial(bins), linewidth=2, color='r')
ax1.set_xlabel("k")
ax1.set_ylabel("P(k)")
ax1.set_title("Histogram of Poisson Distribution")
ax2 = plt.subplot(2, 2, 2)
n, p = 20, 0.5
s = np.random.binomial(n, p, 100000)
count, bins, ignored = ax2.hist(s, 16, normed=True)
ax2.plot(bins, scipy.special.comb(n, bins) * p**bins * (1-p)**(n-bins), linewidth=2, color='r')
ax2.set_xlabel("k")
ax2.set_ylabel("P(k)")
ax2.set_title("Binomial, n=20, p=0.5")
ax3 = plt.subplot(2, 2, 3)
# mu, sigma = 0, 0.1
# s = np.random.normal(mu, sigma, 10000)
# count, bins, ignored = ax3.hist(s, 30, density=True)
mu, sigma = 0, 0.1 # mean and standard deviation
s = np.random.normal(mu, sigma, 1000)
count, bins, ignored = ax3.hist(s, 30, density=True)
ax3.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) *
np.exp( - (bins - mu)**2 / (2 * sigma**2) ),
linewidth=2, color='r')
ax3.set_xlabel("x")
ax3.set_ylabel("P(x)")
ax3.set_title("Normal, mu={}, sigma={}".format(mu, sigma))
ax4 = plt.subplot(2, 2, 4)
a, m = 3, 2. # shape and mode
s = (np.random.pareto(a, 1000) + 1) * m
count, bins, _ = ax4.hist(s, 100, density=True)
fit = a*m**a / bins**(a+1)
ax4.plot(bins, max(count)*fit/max(fit), linewidth=2, color='r')
ax4.set_xlabel("x")
ax4.set_ylabel("P(x)")
ax4.set_title("Pareto")
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment