Scribd
Upload
8 views12 pages

Probability Distributions and Analysis

The document contains Python code for statistical analysis using libraries such as matplotlib, numpy, and scipy. It includes simulations of random outcomes, calculations of probabilities for odd outcomes, and visualizations of various probability distributions including gamma, exponential, Poisson, and normal distributions. Additionally, it analyzes Titanic passenger data to estimate and visualize the probability mass functions for survival, class, gender, and age.

Uploaded by

nishtharakan
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
8 views12 pages

Probability Distributions and Analysis

The document contains Python code for statistical analysis using libraries such as matplotlib, numpy, and scipy. It includes simulations of random outcomes, calculations of probabilities for odd outcomes, and visualizations of various probability distributions including gamma, exponential, Poisson, and normal distributions. Additionally, it analyzes Titanic passenger data to estimate and visualize the probability mass functions for survival, class, gender, and age.

Uploaded by

nishtharakan
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

3/2/25, 12:17 AM project

In [34]: import [Link] as plt


import numpy as np
import random
import [Link] as stats
import seaborn
import pandas as pd

Problem 1.
In [35]: # a.
flips = [10,50,100,500,1000]
x = [1,2,3,4,5]

for i in range(5):

data = [Link](5)
for j in range(flips[i]):
result = [Link](1,5)
data[result-1]+=1
totalodd = float(data[0]+data[2]+data[4])
print("Number of tosses: "+str(flips[i])+". Estimated probabilty of odd

Number of tosses: 10. Estimated probabilty of odd outcome 0.6


Number of tosses: 50. Estimated probabilty of odd outcome 0.5
Number of tosses: 100. Estimated probabilty of odd outcome 0.56
Number of tosses: 500. Estimated probabilty of odd outcome 0.6
Number of tosses: 1000. Estimated probabilty of odd outcome 0.603

1B. There are five possible outcomes, and three are odd (desired), so P(odd)=3/5
1c. Yes, the distribution matches what I would expect, especially for when there are more
rolls.
1D. Next cell
In [36]: #1d
# Chance of 3,4, and 5 are equal, and each is 1/2 of 1 and 2.
# with this info, calculate chance of 3/4/5 as 1/7, and chance of 1 and 2 as
# Predicted ratio of odd rolls = 2/7 + 1/7 + 1/7 ≈ 0.5714

# Once again for higher # of rolls the results are a closer match to the the

for i in range(5):
data = [Link](5)
for j in range(flips[i]):
result = [Link](x,weights=(2,2,1,1,1),k=1)
data[result[0]-1]+=1

# axs[i].bar(x,data)

[Link] 1/12
3/2/25, 12:17 AM project

# axs[i].set_xticks(x)
# fig.set_figheight(15)
totalodd = float(data[0]+data[2]+data[4])
print("Number of tosses: "+str(flips[i])+". Estimated probabilty of odd

Number of tosses: 10. Estimated probabilty of odd outcome 0.5


Number of tosses: 50. Estimated probabilty of odd outcome 0.46
Number of tosses: 100. Estimated probabilty of odd outcome 0.59
Number of tosses: 500. Estimated probabilty of odd outcome 0.544
Number of tosses: 1000. Estimated probabilty of odd outcome 0.566

Problem 2
In [37]: alphas= [0.5,1,1.5,2.5]
lambdaa = 0.5
x = [Link](start=0, stop=10, step=0.1) # 0 to 10
gamma = [Link]

fig,ax = [Link](2)

line1, = ax[0].plot(x,[Link](x=x,a=alphas[0],scale = 1/lambdaa),color = '


#ax[0].plot([Link](x=x,a=alphas[0],scale = 1/lambdaa),color = 'purple')
line2, =ax[0].plot(x,[Link](x=x,a=alphas[1],scale = 1/lambdaa),color = 'r
line3, =ax[0].plot(x,[Link](x=x,a=alphas[2],scale = 1/lambdaa),color = 'g
line4, =ax[0].plot(x,[Link](x=x,a=alphas[3],scale = 1/lambdaa),color = 'g

#ax[0].set_xlim(0,10)
ax[0].set_title( "gamma CDF")
ax[0].set_ylabel("F(x)")
ax[0].set_xlabel("x")
#[Link]()

ax[1].plot(x,[Link](x=x,a=alphas[0],scale = 1/lambdaa),color = 'blue')


#ax[1].plot([Link](x=x,a=alphas[0],scale = 1/lambdaa),color = '')
ax[1].plot(x,[Link](x=x,a=alphas[1],scale = 1/lambdaa),color = 'red')
ax[1].plot(x,[Link](x=x,a=alphas[2],scale = 1/lambdaa),color = 'green')
ax[1].plot(x,[Link](x=x,a=alphas[3],scale = 1/lambdaa),color = 'gray')
#ax[1].set_xlim(0,10)
ax[1].set_title( "gamma PDF")

ax[1].set_ylabel("f(x)")
ax[1].set_xlabel("x")

[Link]([line1,line2,line3,line4],["a=0.5","a=1.0","a=1.5","a=2.5"])
fig.tight_layout()
ax[0].grid()
ax[1].grid()
#[Link](x,loc=0,scale=1)
#[Link](x,loc=0,scale=2)

[Link] 2/12
3/2/25, 12:17 AM project

In [38]: #2b. Exponential

fig,ax = [Link](2)
x = [Link](start=0, stop=10, step=0.1) # 0 to 10
exp = [Link]

ax[0].plot(x,[Link](x=x,scale=1))
ax[0].set_title("Exponential CDF")
ax[0].set_ylabel("F(x)")
ax[0].set_xlabel("x")
#ax[0].set_xlim(0,10)

ax[1].plot(x,[Link](x=x,scale=1))
ax[1].set_title("Exponential PDF")
#ax[1].set_xlim(0,10)
ax[0].grid()

ax[1].set_ylabel("f(x)")
ax[1].set_xlabel("x")
ax[1].grid()
fig.tight_layout()

[Link] 3/12
3/2/25, 12:17 AM project

In [39]: #Poisson
fig,ax = [Link](2)
x = [Link](start=-1, stop=10, step=0.1) # 0 to 10
poisson = [Link]
ax[0].plot(x,[Link](x,mu=3))

ax[0].set_title("Poisson CDF")
ax[0].set_ylabel("F(x)")
ax[0].set_xlabel("x")

ax[1].plot([Link](0,10),[Link]([Link](0,10),mu=3),linestyle = "No
ax[1].set_title("Poisson PMF") # PMF becuase
ax[1].set_ylabel("P(x)")
ax[1].set_xlabel("x")
ax[0].grid()
ax[1].grid()
fig.tight_layout()

[Link] 4/12
3/2/25, 12:17 AM project

In [40]: #2b. Normal Distribution

fig,ax = [Link](2)
x2 = [Link](-5,8,0.1)
norm = [Link]
line1, =ax[0].plot(x2,[Link](x=x2,scale=1))
line2, =ax[0].plot(x2,[Link](x=x2,loc=3,scale=[Link](2)))
ax[0].set_title("Normal Distribution CDF")
ax[0].set_ylabel("F(x)")
ax[0].set_xlabel("x")
#ax[0].set_xticklabels([Link](-10, 10, 1))

#ax[0].set_xlim(0,10)

ax[1].plot(x2,[Link](x=x2))
ax[1].plot(x2,[Link](x=x2,loc=3,scale=[Link](2)))
ax[1].set_title("Normal Distribution PDF")
ax[1].set_ylabel("f(x)")
ax[1].set_xlabel("x")

ax[0].grid()
ax[1].grid()

ax[0].legend([line1,line2],["Standard", "Modified"])
ax[1].legend([line1,line2],["Standard", "Modified"])

[Link] 5/12
3/2/25, 12:17 AM project

fig.set_figheight(10)

fig.tight_layout()

[Link] 6/12
3/2/25, 12:17 AM project

[Link] 7/12
3/2/25, 12:17 AM project

Problem 3 The price class, gender (1 for Male and 0 for Female) and age are recorded
for each sur- vived (1) or killed (0) passenger. We use random variables S, C, G, and A
for survival status, price class, gender and age, respectively
In [41]: #(a) Estimate the PMFs for S, C, G, and A by finding the fraction of each re
#of these random variables among all data. Plot these PMFs.
titanic = pd.read_excel("modified_titanic.xlsx")

fig,ax = [Link](4)

total = titanic["Age"].count()

ax[0].bar(x=[0,1],height=[titanic["Survived"].value_counts().get(0)/titanic[
ax[0].set_ylim(0,1)
ax[0].set_ylabel("P(S)")
ax[0].set_xlabel("Survival")
ax[0].set_title("PMF of survived")
ax[0].set_xticks(ticks=[0,1],labels=["Dead","Survived"])

ax[1].bar(x=[1,2,3],height=[titanic["Pclass"].value_counts().get(1)/titanic[
ax[1].set_ylim(0,1)
ax[1].set_ylabel("P(C)")
ax[1].set_xlabel("Class")
ax[1].set_title("PMF of class")
ax[1].set_xticks(ticks=[1,2,3])

ax[2].bar(x=[0,1],height=[titanic["Sex"].value_counts().get(0)/titanic["Sex"
ax[2].set_ylim(0,1)
ax[2].set_ylabel("P(G)")
ax[2].set_xlabel("Gender")
ax[2].set_title("PMF of gender")
ax[2].set_xticks(ticks=[0,1],labels=["Female","Male"])

age_data = [Link](90)
#print(titanic["Age"].value_counts().get(0))
for i in range(90):
if (titanic["Age"].value_counts().get(i)!=None):
age_data[i] = titanic["Age"].value_counts().get(i)/titanic["Age"].co
ax[3].bar(x=[Link](0,90),height = age_data)
ax[3].set_title("PMF of age")
ax[3].set_ylabel("P(A)")
ax[3].set_xlabel("Age")

y_pos =[Link](len(titanic["Age"].unique()))
fig.set_figheight(8)
fig.tight_layout()

[Link] 8/12
3/2/25, 12:17 AM project

In [42]: #3b
fig,ax = [Link](3,layout = 'constrained')
alive = titanic["Survived"] == 1
dead = titanic["Survived"] == 0
totaldead = [Link]()
totalalive=[Link]()

# Code is from documentation on this type of graph


#[Link]
conditional_sex = {

[Link] 9/12
3/2/25, 12:17 AM project

'Dead':[(dead &(titanic["Sex"]==0)).sum()/totaldead,(dead&(titanic["Sex"
"Survived":[(alive &(titanic["Sex"]==0)).sum()/totalalive,(alive&(titani
}
x = [Link](2)
width = 0.25
multiplier = 0

for attribute, measurement in conditional_sex.items():


offset = width*multiplier
rects = ax[0].bar(x + offset, measurement,width, label=attribute)
ax[0].bar_label(rects, padding=3)
multiplier += 1
ax[0].set_ylabel('P(G|S)')
ax[0].set_title('Sex conditioned by Survival')
ax[0].set_xticks(ticks=[0+width/2,1+width/2],labels=["Female","Male"])
ax[0].legend()
ax[0].set_ylim(0,1)

conditional_class={'Dead':[],"Survived":[]}
for i in range(3):
conditional_class["Dead"].append((dead &(titanic["Pclass"]==i+1)).sum()/
conditional_class["Survived"].append((alive &(titanic["Pclass"]==i+1)).s

x = [Link](3)
width = 0.25
multiplier = 0

for attribute, measurement in conditional_class.items():


offset = width*multiplier
rects = ax[1].bar(x + offset, measurement,width, label=attribute)
ax[1].bar_label(rects, padding=3)
multiplier += 1
ax[1].set_ylabel('P(C|S)')
ax[1].set_title('Class conditioned by Survival')
ax[1].set_xticks(ticks=x+width/2,labels=["1","2",'3'])
ax[1].legend()
ax[1].set_ylim(0,1)

x = [Link](90)
conditional_age={'Dead':[],"Survived":[]}
for i in x:
conditional_age["Dead"].append((dead &(titanic["Age"]==i)).sum()/totalde
conditional_age["Survived"].append((alive &(titanic["Age"]==i)).sum()/to

width = 0.25
multiplier = 0

for attribute, measurement in conditional_age.items():


offset = width*multiplier
rects = ax[2].bar(x + offset, measurement,width, label=attribute)
#ax[2].bar_label(rects, padding=3)
multiplier += 1
ax[2].set_ylabel('P(A|S)')
ax[2].set_title('Age conditioned by Survival')
ax[2].set_xticks(ticks=[Link](0,90,10))
ax[2].legend()

[Link] 10/12
3/2/25, 12:17 AM project

Out[42]: <[Link] at 0x169b0ebd0>

In [43]: #3c
#we first do the one where S =0
p_age_dead=0
for i in range(41): # range from 0 to 40
p_age_dead+=conditional_age["Dead"][i]
p_class_dead = conditional_class["Dead"][1 - 1]
p_sex_dead = conditional_sex["Dead"][0]
p_AB_dead = p_age_dead*p_class_dead*p_sex_dead * totaldead/total
print("P(S = 0,C = 1,G= 0,A≤40) = " + str(p_AB_dead)) #last term is P(S=0)

p_age_alive=0
for i in range(41): # range from 0 to 40
p_age_alive+=conditional_age["Survived"][i]
p_class_alive = conditional_class["Survived"][0]
p_sex_alive = conditional_sex["Survived"][0]
p_AB_alive = p_age_alive*p_class_alive*p_sex_alive * totalalive/total
print("P(S = 1,C = 1,G= 0,A≤40) = " + str(p_AB_alive)) #last term is P(S=1)

P(S = 0,C = 1,G= 0,A≤40) = 0.0103547704894279


P(S = 1,C = 1,G= 0,A≤40) = 0.08216205203987068

In [44]: #3d
#P(A|B) = P(B|A)*P(A)/P(B), so P(B|A)=P(A|B)*P(B)/P(A) = P(A and B) / P(A),
# A here is the combination of C G A, and B is S.
print("P(S = 0|C = 1,G= 0,A≤40) death = " + str(p_AB_dead /(p_age_dead*p_cla
print("P(S = 1|C = 1,G= 0,A≤40) survival = " + str(p_AB_alive /(p_age_alive*

[Link] 11/12
3/2/25, 12:17 AM project

P(S = 0|C = 1,G= 0,A≤40) death = 0.6144306651634723


P(S = 1|C = 1,G= 0,A≤40) survival = 0.38556933483652767

The probability of dying is greater than the probability of surviving, so I predict dying.
In [45]: #4A

test_n = [1,3,10,30,100]

for k in test_n:
n=k
z = 0
for i in range(n):
x = [Link]()*4+3
z+=x
z/=n

[Link] 12/12

You might also like