iframe-proxy | Sunbelt Computer Software

155 lines (117 loc) · 4.23 KB
Created on Mon Jun 26 00:55:07 2019
@author: Neeraj
Description: This code illustrates how to do basic porbability operations in Python. 
Reference: Chapter 6 : Probability
import enum, random
class Kid(enum.Enum):
def random_kid() -> Kid:
  return random.choice([Kid.BOY, Kid.GIRL])
both_girls = 0
older_girl = 0
either_girl = 0
random.seed(0) # Ensures reproducibility of the experiments
for _ in range(10000):
  younger = random_kid()
  older = random_kid()
  if older == Kid.GIRL:
    older_girl += 1
  if older == Kid.GIRL and younger == Kid.GIRL:
    both_girls += 1
  if older == Kid.GIRL or younger == Kid.GIRL:
    either_girl += 1
print(f"P(both girls | older girl = {both_girls/older_girl})")
print(f"P(both girls | either girl = {both_girls/either_girl})")
# Create pdf of a uniform random variable
def uniform_pdf(x: float) -> float:
  """A uniform random variable gives equal probability to values between 0 and 1."""
  return 1 if  0 <= x < 1 else 0
# Create cdf of a uniform random variable
def uniform_cdf(x: float) -> float:
  if x <= 0: return 0 
  elif x < 1: return x
  else: return 1
# Plot uniform pdf
import numpy as np
x = np.linspace(-0.5,1.5,21)
for i in x:
    y.append(uniform_pdf(i))
import matplotlib.pyplot as plt
plt.plot(x,y)
# Plot uniform cdf
for i in x:
    y.append(uniform_cdf(i))
import matplotlib.pyplot as plt
plt.plot(x,y)
plt.axis([-0.5, 1.5, -0.5, 1.2])
# The normal distribution
from math import sqrt, pi, exp
def normal_pdf(x: float, mu: float = 0, sigma: float = 1) -> float:
    num = exp(-(x-mu)**2 / 2 / sigma**2)
    den = sqrt(2*pi)*sigma
    return num/den
# Plot normal pdf
x = [x/10.0 for x in range(-50,50)]
y = [normal_pdf(i) for i in x]
plt.plot(x,y, '-', label ='mu = 0, sigma = 1')
plt.plot(x,[normal_pdf(i,0,2) for i in x], '--', label ='mu = 0, sigma = 2')
plt.plot(x,[normal_pdf(i,1,1) for i in x], '-.', label ='mu = 1, sigma = 1')
from math import erf
def normal_cdf(x:float, mu: float = 0, sigma: float = 1) -> float:
    return (1 + erf((x-mu)/sqrt(2) / sigma))/2
x = [x/10.0 for x in range(-50,50)]
plt.plot(x,[normal_cdf(i,0,1) for i in x], '-', label ='mu = 0, sigma = 1')
plt.plot(x,[normal_cdf(i,0,2) for i in x], '--', label ='mu = 0, sigma = 2')
plt.plot(x,[normal_cdf(i,1,1) for i in x], '-.', label ='mu = 1, sigma = 1')
# Inverse of normal distribution
def inverse_normal_cdf(p: float, 
                       mu: float = 0, 
                       sigma: float = 1,
                       tol: float = 0.01) -> float:
    """Computes a value from Z(0,1) at specified probability (p) level"""
    # Convert to standard normal distribution
    if mu != 0 or sigma != 1:
        return mu + sigma * inverse_normal_cdf(p)
    low_z = -10.0 # standard normal is 0 at cdf(-10)
    high_z = 10.0 # standard normal is 0 at cdf(10)
    # Use binary search to find the desired value
    while high_z - low_z > tol:
        midz = (high_z + low_z)/2
        pmid = normal_cdf(midz)
        if p > pmid:
            low_z = midz
        else:
            high_z = midz
        #print(midz)
    return midz
print(inverse_normal_cdf(0.3))      
import random
def bernoulli_trial(p: float) -> int:
    """Returns 1 with probability p and 0 with probability 1-p"""
    return 1 if random.random() < p else 0
def binomial(n: int, p: float) -> int:
    """Returns the sum of n bernoulli(p) trials"""
    return sum(bernoulli_trial(p) for _ in range(n))
from collections import Counter
def binomial_histogram(p: float, n: int, num_points: int) -> None:
    """Picks points from a Binomial(n, p) and plots their histogram"""
    data = [binomial(n, p) for _ in range(num_points)]
    # use a bar chart to show the actual binomial samples
    histogram = Counter(data)
    plt.bar([x - 0.4 for x in histogram.keys()],
            [v / num_points for v in histogram.values()],
            color='0.75')
    mu = p * n
    sigma = sqrt(n * p * (1 - p))
    # use a line chart to show the normal approximation
    xs = range(min(data), max(data) + 1)
    ys = [normal_cdf(i + 0.5, mu, sigma) - normal_cdf(i - 0.5, mu, sigma)
          for i in xs]
    plt.plot(xs,ys)
    plt.title("Binomial Distribution vs. Normal Approximation")
Sunbelt Computer Software

PL/B Language Development and Support

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Sunbelt Computer Software

PL/B Language Development and Support

FilesExpand file tree

probability.py

Latest commit

History

probability.py

File metadata and controls
