<jupyter_start><jupyter_text>Modeling Distributions and RNGsSo, last week we learned about PMFs, PDFs, CDFs and all that. This week we're going to learn how to *use* the relationship between CDFs and PDFs to transform a "random number generator" (RNG) from one distribution to another. It's cool! Let's review the general properties of a CDF to get things started.A CDF for a random variable $t$ is defined as:$$\text{CDF(t) = Probability (random var is between $-\infty$ and $t$)} = \int_{-\infty}^{t} P(t')\,dt'$$In other words it's the probability that the random variable is measured with a value less than $t$.The domain of the CFD is $-\infty$ to $\infty$ and the range is between 0 and 1, since the value of the function is a probability. This means that *any* possible value of the random variable can be used as an argument, and the result will *always* be a probability between 0 and 1. Let's take a concrete example, like the one we used last week:$$ PDF(t) = \frac{1}{\tau} e^{-\frac{t}{\tau}} $$$$ CDF(t) = 1-e^{-\frac{t}{\tau}} $$Let's look at the graph of the CDF:<jupyter_code>%matplotlib inline
import matplotlib.pyplot as pl
import matplotlib as mp
import numpy as np

t=np.linspace(0,70,200)
Tau=10.0

def CDF(t):
    return 1-np.exp(-t/Tau)

pl.figure(figsize=(10,6))
pl.plot(t, CDF(t), 'b-')
pl.title("CDF(t) for exponential distribution")
pl.xlabel("time (min)")
pl.ylabel("CDF(t)")
pl.grid()<jupyter_output><empty_output><jupyter_text>Now suppose we sprinkle fairy dust uniformly along the vertical axis of this graph. For each speck of dust we can reverse map to the time axis like so:<jupyter_code>N=300 # how many specks of dust?

#
# set up some functions to compute the CDF and inverse CDF
# Don't worry too much about this for the moment, but just focus on the picture.
# Later you can come back and figure out how it works if you like.
#

fdusts = np.random.rand(N) # these are our fairy dust specs
rev = lambda r: -Tau*np.log(1.0-r)
cdf = lambda t: 1.0 - np.exp(-t/Tau)

pl.figure(figsize=(10,6))
pl.plot(t, CDF(t), 'b-')
pl.title("CDF(t) for exponential distribution with fairy dust!")
pl.xlabel("time (min)")
pl.ylabel("CDF(t)")
pl.grid()
currentAxis = pl.gca()
for yval in fdusts:
    xval=rev(yval)
    pl.plot([0,xval],[yval,yval],'g:')
    pl.plot([xval,xval],[yval,0],'g:')<jupyter_output><empty_output><jupyter_text>Note that if we move the bits of fairy dust to the time axis by reverse mapping the y values from the CDF backwards to the corresponding times. For any given chunk of CDF range there there is a proporational number of bits of fairy dust on the time axis. What happens to the *density* of fairy dust on the "time" axis? Where is it most dense? Least? If you think about it for a bit, you'll see that the fairy dust piles up on the time axis with a density that is proportional to the *slope* of the CDF. For any given "$\Delta t$" the amount of "$\Delta P$" (fairy dust) is equal to:$$ \Delta P = \left(\frac{d(CDF(t))}{dt}\right) \Delta t $$But what is this derivative? Think about it. ;-) Yes, it's nothing other than the $PDF(t)$!So by following this recipe, (back mapping a random number between 0 and 1 through the CDF) we get a distribution of *times* that is distributed exactly as the PDF requires. Simple Recipe for RNGs1) Start with the PDF, and find the CDF$$\text{CDF(t) = Probability (random var is between $-\infty$ and $t$)} = \int_{-\infty}^{t} P(t')\,dt'$$2) Solve the following equation for $t$:$$ r = CDF(t);\ \ \ \ \ \ (Eq 1)$$where $r$ is a uniformly distributed random number between 0 and 1.3) Use the values of $t$ as the random variable.Of course, this only works if you can actually *solve* for $t$ in ${\rm Eq.}\ 1$. An example: The Exponential DistributionLet's try this with the exponentail distribution where we already know the CDF:$$ CDF(t) = 1.0 - e^{-t/\tau} $$For step 2:$$ r = 1.0 - e^{-t/\tau}$$so:$$ e^{-t/\tau} = 1.0 - r $$$$ -\frac{t}{\tau} = \ln(1.0-r) $$$$ t = - \tau \ln(1.0-r) $$OK! So we can generate a bunch ($N$) of 'r' numbers and transform them into a bunch of exponentially distributed values of 't' according to the exponential PDf. Let's try it!<jupyter_code>N=10000 # how many random numbers to generate?

r = np.random.rand(N)
t = -Tau*np.log(1.0-r)

n,_,_ =pl.hist(t, bins=np.linspace(0,50,41))
pl.title("Exponentially distributed 't's")
pl.xlabel("t (min)")
pl.ylabel("Number (out of %d)" % N)

pl.plot(np.linspace(0,50,41)[:-1],np.log(n),'b.')
pl.grid()<jupyter_output><empty_output><jupyter_text>Once you have this distribution you can begin to ask other kinds of interesting questions. For example, suppose we wanted to estimate the chances of seeing 4 or more defects in less than 10 minutes. Or, we need to estimate how big of a backlog we're likely to see by the end of the day if our defect rate when up by 10% for 3 hours? These are questions that might be tricky to handle analytically, but can easily be simulated using RNGs.There are also some python "tricks" that you might find useful in doing these kinds of simulations.Suppose we simulate a set of intervals:<jupyter_code>t = -Tau*np.log(1.0-np.random.rand(100))  # get 100 intervals between defects
pl.hist(t)

# What is the total time for these 100 intervals?

print("Total time:", t.sum())
# What is the averate time between defects for these samples?

print("Average time:", t.mean())
# What is the time between four successive defects? Use the "slice" concept

delta3 = t[:-2] + t[1:-1] + t[2:]
pl.hist(delta3)
print("Minimum time between four successive defects:", min(delta3))

ThreeInLessThan10 = delta3<10
print("times less than 10:", delta3[np.nonzero(ThreeInLessThan10)])
print("Number of times 4 defects happened in less than 10 minutes:", len(delta3[np.nonzero(ThreeInLessThan10)]))
#
# Let's simulate an assembly line for one hour as an example.
#

finalTime = 60            # run for one hour
processTime = 10.0        # how long does it take to process an 
currTime = 0.0            # what time is it now?
nextTime = processTime    # when will the last defect be finished? (it takes 10 minutes to handle a defect)
defectsProcessed = 0      # how many defects have been processed?
defectsOccurred = 0

for delta in t:
    currTime = currTime + delta # when does the next defect appear?

    if currTime > finalTime:
        break # OK, it's the end of the 1 hour shift, break out.

    defectsOccurred += 1  # count the number of defects that have happened.

    if currTime > nextTime:
        nextTime = currTime + processTime # process immediately
    else:
        nextTime = nextTime + processTime # process after this one.
        
    if nextTime<finalTime:
        defectsProcessed += 1
        
print("Defects Occurred:", defectsOccurred)
print("Defects Processed:", defectsProcessed)
print("Backlog:", defectsOccurred-defectsProcessed)<jupyter_output>Defects Occurred: 5
Defects Processed: 4
Backlog: 1<jupyter_text>Suppose we collect a bunch of "time between defects" data and we want to analyze it to determine if it approximates an exponetial distribution. Can we use a Bayesian approach to determining the parameters of the distribution? Sure!See the statistical exercises for project 4 for an in-depth intro to pymc3, but this is actually pretty quick.<jupyter_code>import pymc3 as pm
import arviz as az

with pm.Model() as model:
    tau = pm.Uniform('tau',lower=0,upper=100)             # set up a prior for tau, somewhere between 0 and 100
    t_obs = pm.Exponential('t_obs',lam=1/tau, observed=t) # add the observed data
    trace = pm.sample(1000, return_inferencedata=True)    # check the results
az.plot_trace(trace)
az.summary(trace, kind="stats")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11300390117035111, 'fraction_numerical': 0.01846553966189857, 'mean_word_length': 4.193112761647535, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '132109', 'license': 'MIT', 'n_tokens_mistral': 2571, 'n_tokens_neox': 2410, 'n_words': 1044, 'path': 'proj3/Modeling Distributions and RNGs.ipynb', 'repo_name': 'sspickle/instrumentation-projects', 'size': '7690'}