### Statistics 6401, Handout 1

### This Python script simulates flipping a fair coin 1000 times.
### Three examples of using simulation are given using various 
### Python commands.  Finally, the CLT is demonstrated.

# Example 1. Simulate flipping a fair coin n times.

import random

n = 1000
p = 0.50

h = 0 	# counter for heads

for _ in range(n):
  x = random.random()
  if x < p:
	   h += 1

print(h)	# total number of heads simulated

ph = h/n
print(ph)	# simulated estimate of the probability of heads

# Example 2. Simulate flipping a fair coin n times.
# Simulate using vectors.


import numpy as np

n = 1000
p = 0.50

x = np.random.uniform(0,1,n)	# creates a vector x of random uniform values

h = sum(x < p)	# gives the number of heads

ph = h/n

print(h)	# total number of heads simulated

print(ph)	# simulated estimate of the probability of heads

# Example 3. Recall that flipping a fair coin n times is a binomial experiment

n = 1000
p = 0.50

h = np.random.binomial(n,p)	# number of heads

ph = h/n

print(h)	# total number of heads simulated

print(ph)	# simulated estimate of the probability of heads


# Example 4. Now simulate flipping a fair coin n times, but do it 100,000 
# times and then plot plot the sampling distribution of the estimated 
# probability of heads.  This is an example of the Central Limit Theorem.

import matplotlib.pyplot as plt

n = 1000
p = 0.50

Reps = 100000

ph = np.zeros(Reps)

for i in range(Reps):
  	h = np.random.binomial(n,p)
  	ph[i] = h/n
	
	
#plt.clf()	
plt.hist(ph,bins=10)
plt.xlim(.45, .55)
plt.ylim(top=40000)
plt.show()
  
####################################################################

### This is an Python program that simulates a system of 
### n components connected in series.  The probability that
### a component fails is p.  The probability that
### the system fails is approximated using simulation.

import random

p = 0.25
n = 4

Reps = 10000

Count2 = 0	# Counter for number of failing systems in the overall simulation
for j in range(Reps):
	Count1 = 0	# Counter for number of failures in a simulated system
	for i in range(n):
		ci = random.random()
		if ci < p:
			Count1 += 1
	if Count1 > 0:
		Count2 += 1

psf = Count2/Reps

print(psf)

psf_truth = 1-(1-p)**n

print(psf_truth)

### Alternative solution, using Python

c = [random.random() < p for _ in range(Reps*n)]

Count = [sum(c[i:i+n]) for i in range(0,Reps*n,n)]

Count = [x for x in Count if x > 0]

psf = len(Count)/Reps

print(psf)

print(psf_truth)

####################################################################