import pandas as pd
import matplotlib.pyplot as pl
import numpy as np
f = [100, 500, 1000, 2000, 5000, 7000, 8000, 9000, 10000]
amp = [976, 800, 560, 320, 160, 128, 120, 104, 104]
g = [2, 1.639, 1.148, 0.656, 0.328, 0.262, 0.246, 0.213, 0.213]
omega = 2*np.pi*np.array(f)
R=1000
C=.22e-6
X = 1.0/(omega*C)
Z = np.sqrt(R**2 + X**2)
G = 2.0*X/Z
pl.loglog(f, G,'r-',label='theory')
pl.loglog(f,g,'b.',label='exp')
pl.grid()
pl.legend()
G
N=10000
Tau=10.0
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()
t = -Tau*np.log(1.0-np.random.rand(100)) # get 100 intervals between defects
pl.hist(t)
print("Total time:", t.sum())
print("Average time:", t.mean())
Total time: 1119.9999928835432
Average time: 11.199999928835432
delta2 = t[:-1] + t[1:]
pl.hist(delta2)#make it delta 2 and sum of two slices not three - Tariq
print("Minimum time between three successive defects:", min(delta2))
ThreeInLessThan10 = delta2<10#check for delta 2 less then ten - Tariq
print("times less than 10:", delta2[np.nonzero(ThreeInLessThan10)])
print("Number of times 3 defects happened in less than 10 minutes:", len(delta2[np.nonzero(ThreeInLessThan10)]))
Minimum time between three successive defects: 3.3049130530015223
times less than 10: [5.14075085 3.30491305 9.92972089 9.57382404 8.41053916 4.00621959
6.14685499 3.53436547 4.63811215 4.19449225 8.16541285 6.3077787
5.64708998 5.31746605 5.91256128 3.71228103]
Number of times 3 defects happened in less than 10 minutes: 16
finalTime = 480 # 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)
Defects Occurred: 42
Defects Processed: 40
Backlog: 2
!pip install pymc3
Requirement already satisfied: pymc3 in /root/venv/lib/python3.7/site-packages (3.11.2)
Requirement already satisfied: patsy>=0.5.1 in /root/venv/lib/python3.7/site-packages (from pymc3) (0.5.1)
Requirement already satisfied: arviz>=0.11.0 in /root/venv/lib/python3.7/site-packages (from pymc3) (0.11.2)
Requirement already satisfied: dill in /shared-libs/python3.7/py/lib/python3.7/site-packages (from pymc3) (0.3.3)
Requirement already satisfied: numpy>=1.15.0 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from pymc3) (1.19.5)
Requirement already satisfied: fastprogress>=0.2.0 in /root/venv/lib/python3.7/site-packages (from pymc3) (1.0.0)
Requirement already satisfied: pandas>=0.24.0 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from pymc3) (1.2.3)
Requirement already satisfied: semver in /root/venv/lib/python3.7/site-packages (from pymc3) (2.13.0)
Requirement already satisfied: cachetools>=4.2.1 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from pymc3) (4.2.1)
Requirement already satisfied: scipy>=1.2.0 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from pymc3) (1.6.1)
Requirement already satisfied: typing-extensions>=3.7.4 in /shared-libs/python3.7/py-core/lib/python3.7/site-packages (from pymc3) (3.7.4.3)
Requirement already satisfied: theano-pymc==1.1.2 in /root/venv/lib/python3.7/site-packages (from pymc3) (1.1.2)
Requirement already satisfied: six in /shared-libs/python3.7/py-core/lib/python3.7/site-packages (from patsy>=0.5.1->pymc3) (1.15.0)
Requirement already satisfied: matplotlib>=3.0 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from arviz>=0.11.0->pymc3) (3.3.4)
Requirement already satisfied: netcdf4 in /root/venv/lib/python3.7/site-packages (from arviz>=0.11.0->pymc3) (1.5.6)
Requirement already satisfied: setuptools>=38.4 in /root/venv/lib/python3.7/site-packages (from arviz>=0.11.0->pymc3) (47.1.0)
Requirement already satisfied: xarray>=0.16.1 in /root/venv/lib/python3.7/site-packages (from arviz>=0.11.0->pymc3) (0.17.0)
Requirement already satisfied: packaging in /shared-libs/python3.7/py-core/lib/python3.7/site-packages (from arviz>=0.11.0->pymc3) (20.9)
Requirement already satisfied: pytz>=2017.3 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from pandas>=0.24.0->pymc3) (2021.1)
Requirement already satisfied: python-dateutil>=2.7.3 in /shared-libs/python3.7/py-core/lib/python3.7/site-packages (from pandas>=0.24.0->pymc3) (2.8.1)
Requirement already satisfied: filelock in /root/venv/lib/python3.7/site-packages (from theano-pymc==1.1.2->pymc3) (3.0.12)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /shared-libs/python3.7/py-core/lib/python3.7/site-packages (from matplotlib>=3.0->arviz>=0.11.0->pymc3) (2.4.7)
Requirement already satisfied: pillow>=6.2.0 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from matplotlib>=3.0->arviz>=0.11.0->pymc3) (8.1.1)
Requirement already satisfied: cycler>=0.10 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from matplotlib>=3.0->arviz>=0.11.0->pymc3) (0.10.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from matplotlib>=3.0->arviz>=0.11.0->pymc3) (1.3.1)
Requirement already satisfied: cftime in /root/venv/lib/python3.7/site-packages (from netcdf4->arviz>=0.11.0->pymc3) (1.4.1)
WARNING: You are using pip version 20.1.1; however, version 21.0.1 is available.
You should consider upgrading via the '/root/venv/bin/python -m pip install --upgrade pip' command.
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
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [tau]
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 3 seconds.
az.plot_trace(trace)
az.summary(trace, kind="stats")