AMS326HW3Q2Spring2020.py

#Ryan Burgert AMS326 HW3 Q2 Spring 2020
#   I used 1984444 trials instead of 19844444444444 like specified in the homework, because it gives the same answers (albeit at slightly less precision)
#   Originally this code was meant to calculate the probability of needles being thrown onto the table; I modified it once I realized the question was asking for discs (and not needles)
#   If you want to see how it works with needles, replace and relpace this code, replacing "type=disc" with "type=needle"
#   SAMPLE OUTPUT SHOWN BELOW:
#      Below is a table that shows the proability of hitting exactly k lines given a disc of diameter d and line spacing w=1
#             d=0.75     d=0.1     d=0.3     d=0.5     d=0.6     d=0.7     d=0.8     d=0.9     d=1.5  d=2.0  d=3.0
#      k=0  0.249757  0.900294  0.699491  0.500899  0.399538  0.300020  0.200081  0.100059  0.000000    0.0    0.0
#      k=1  0.749574  0.100025  0.299628  0.499491  0.600155  0.699893  0.800064  0.900145  0.500391    0.0    0.0
#      k=2  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.499718    1.0    0.0
#      k=3  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000    0.0    1.0
#      k=4  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000    0.0    0.0
#A link to a gist of this code: 

from rp import *#pip install rp

def random_theta_values(quantity=10000):
    from numpy.random import rand
    return rand(quantity)*2*pi
def zone(x,line_width):
    from numpy import floor
    return floor(x/line_width)
def random_x_positions(line_width=1,quantity=10000):
    from numpy.random import rand
    return rand(quantity)*line_width
def how_many_lines_does_needle_cross(needle_x,line_width,theta,needle_length):
    from numpy import cos,abs
    squashed_needle_width=abs(cos(theta)*needle_length)
    needle_left_x =needle_x-squashed_needle_width/2
    needle_right_x=needle_x+squashed_needle_width/2
    needle_left_zone =zone(needle_left_x ,line_width)
    needle_right_zone=zone(needle_right_x,line_width)
    return np.abs(needle_right_zone-needle_left_zone)#the needle belongs to two zones
def estimate_probability(quantity=10000,line_width=1,needle_length=1,k=1,type='disc'):
    #Probability that we cross exactly k lines
    assert type in ['disc','needle']
    needle_xs      =random_x_positions    (line_width,quantity                      )
    thetas         =random_theta_values   (quantity                                 )
    if type=='disc':
        #A disc with diameter d is equivalent to a horizontal needle with length d
        thetas*=0#We're using discs, which is equivalent to a bunch of needles
    hits_and_misses=k==how_many_lines_does_needle_cross(needle_xs,line_width,thetas,needle_length)
    return np.mean(hits_and_misses)

big_number=1984444
def run_test(d=3/4,k=1,w=1):
    p=estimate_probability(quantity=big_number,needle_length=d,line_width=1,k=k)
    print("Probability of crossing exactly",k,'lines with disc diameter',d,'and line spacing',w,'is',p)
    return p

diameters=[3/4,1/10,3/10,5/10,6/10,7/10,8/10,9/10,15/10,20/10,30/10]
k_values=[0,1,2,3,4]


import pandas as pd

df={}
for d in [3/4,1/10,3/10,5/10,6/10,7/10,8/10,9/10,15/10,20/10,30/10]:
    w=1
    row=[]
    for k in [0,1,2,3,4]:
        row.append(run_test(d,k,w))
    df['d='+str(d)]=row

print()
print("Below is a table that shows the proability of hitting exactly k lines given a disc of diameter d and line spacing w=1")
print()

df = pd.DataFrame(df,index=['k='+str(k) for k in k_values])
print(df.to_string())#print a chart in the console