Chapter 5: Probability Densities

Example, page-122

In [1]:
#Variable Declaration

#Calculation
from scipy import *
from pylab import *
from scipy.integrate import *

# between 1 and 3
result = quad(lambda x: 2*pow(e,-2*x), 1, 3)         # probability for 1 < x < 3
result=round(result[0],3)

# greater than 0.5
result1 = quad(lambda x: 2*pow(e,-2*x), 0.5, +inf)   # probability for 0.5 < x < +inf
result1 = round(result1[0],3)

# Result
print "probability of x to lie between 1 and 3 is: ",result
print "probability of x to be greater than 0.5 is: ",result1

probability of x to lie between 1 and 3 is:  0.133
probability of x to be greater than 0.5 is:  0.368

Example, page-123

In [2]:
# Calculation
from scipy import *
from pylab import *
from scipy.integrate import *

# x<=1
result = array(quad(lambda x: 0, -inf, 0)) + array(quad(lambda x: 2*pow(e,-2*x), 0, 1))    # probability for -inf < x <= 1
result = round(result[0],3)

# Result
print "probability of x less than or equal to 1 is: ",result

probability of x less than or equal to 1 is:  0.865

Example, Page-124

In [3]:
#Variable Declaration

#Calculation
from scipy import *
from pylab import *
from scipy.integrate import *

# between 1 and 3
result = quad(lambda x: x*2*pow(e,-2*x), 0, +inf)         # probability for 1 < x < 3
result=round(result[0],3)

# greater than 0.5
result1 = quad(lambda x: pow((x-result),2)*2*pow(e,-2*x), 0, +inf)   # probability for 0.5 < x < +inf
result1 = round(result1[0],3)

# Result
print "Mean: ",result
print "Variance: ",result1

Mean:  0.5
Variance:  0.25

Example, page-126

In [4]:
# Variable Declaration 

# Calculation
from scipy import *
from pylab import *

# (A) between 0.87 and 1.28
# p1 = p( at 1.28 ) - p( at 0.87 )
p1 = 0.8997 - 0.8078              # probability between 0.87 and 1.28
p1 = round(p1,4)

# (B) between -0.34 and 0.62
# p2 = p( at 0.62 ) - p( at -0.34 )
p2 = 0.7324 - 0.3669              # probability between -0.34 and 0.62
p2 = round(p2,4)
 
# (C) greater than 0.85
# p3 = 1 - p( at 0.85 )
p3 = 1 - 0.8023                   # probability greater than 0.85
p3 = round(p3,4)
 
# (D) greater than -0.65
# p4 = 1 - p( at -0.65 )
p4 = 1 - 0.2578                   # probability greater than -0.65
p4 = round(p4,4)

# Result
print "probability 0.87<x1.28 : ",p1
print "probability -0.34<x<0.62 : ",p2
print "probability x>0.85 : ",p3
print "probability x> (-0.65) : ",p4

probability 0.87<x1.28 :  0.0919
probability -0.34<x<0.62 :  0.3655
probability x>0.85 :  0.1977
probability x> (-0.65) :  0.7422

Example, Page-128

In [5]:
# Variable declaration

# Calculation
from scipy import *
from pylab import *

# as we know F(z:0.01) = 0.99 & closest to the value 0.99 i.e. 0.9901 so z = 2.33
# as we know F(z:0.05) = 0.95 & closest to the value 0.95 i.e. 0.9495 & 0.9505 so z = (1.64+1.65)/2

# Result
print "Part(a): ",2.33
print "Part(b): ",(1.64+1.65)/2

Part(a):  2.33
Part(b):  1.645

Example, Page-129

In [6]:
# Variable declaration
Mean = 10.1      # dB
std_dev = 2.7    # dB

# Calculation
from scipy import *
from pylab import *

# Part(a)
Lower = (8.5-Mean)/std_dev 
Upper = (13.0-Mean)/std_dev

# We need to calculate F- Value for lower & Upper values.
# F(1.07) = 0.8577 & F(-0.59) = 0.2776

# Part(b)
val = 1 - 0.9678   # 1 - F(1.85)
# Result
print "Part(a): ",(0.8577-0.2776)
print "Part(b): ",(0.8577-0.2776)
print "Part(c): ",(1-0.9678)

Part(a):  0.5801
Part(b):  0.5801
Part(c):  0.0322

Example, page-130

In [7]:
# Variable declaration
std_dev = 0.04      #standard deviation
x = 4               # number of ounces
p = 0.02            # probability of x <4

# Calculation
from scipy import *
from pylab import *

# as we know p = (x-Mean)/std_dev
# from Table-3 (page-514) we find closest value to 0.02 which is 0.0202 corresponding to Z = -2.05
Mean = 4 - (-2.05*0.04)            # from z=(x-mean)/std_dev
Mean=round(Mean,3)

# Result
print "Mean of given distribution: ",Mean,"ounces"

Mean of given distribution:  4.082 ounces

Example, Page-130

In [8]:
# Variable declaration
Mean = -4.6
variance = 1.21 
x = 0.0015   # oz/st gold

# Calculation
from scipy import *
from pylab import *

Lower = (math.log(x,e)-Mean)/sqrt(variance)

# We need to calculate F- Value for lower i.e. F(-1.73) to be find out.
# Using table F(-1.73) = 0.0419
# As F-value is too small 

# Result
print "Specimen was collected outside of vein."

Specimen was collected outside of vein.

Example, Page-131

In [9]:
# Variable declaration
Mean = 11.6
std_dev = 3.3 

# Calculation
from scipy import *
from pylab import *

Lower = 7.5   # corresponding to outages = 8 -> 7.5 to 8.5 thus atleast = 7.5

index = (Lower-Mean)/std_dev
F = 0.1075     # F value for F(index) i.e. F(-1.24) = 0.1075

# Result
print "Desired probability: ",1-F

Desired probability:  0.8925

Example, Page-132

In [10]:
# Variable declaration
p = 0.2
n = 100 

# Calculation
from scipy import *
from pylab import *

Mean = n*p
std_dev = sqrt(n*p*(1-p))

Upper = 15.5   # corresponding to chips = 15 -> 14.5 to 15.5 thus atmost = 15.5
Lower = 14.5

in_upper = (Upper-Mean)/std_dev
in_lower = (Lower-Mean)/std_dev

f_upper = 0.1292     # F value for F(upper) i.e. F(-1.13) = 0.1292
f_lower = 0.0838     # F value for F(upper) i.e. F(-1.38) = 0.0838

# Result
print "Part(A) Probability: ",f_upper
print "Part(B) Probability: ",(f_upper-f_lower)

Part(A) Probability:  0.1292
Part(B) Probability:  0.0454

Example, page-137

In [11]:
# Variable declaration
alpha = 2               # mean of normal distribution
beta = 0.1              # standard deviation of normal distribution

# Calculation
from scipy import *
from pylab import *

# we need to find p(Io/Ii) = p( (ln(b)-alpha) / beta) - p( (ln(a)-alpha) / beta)
a = 6.1                                        # lower limit
b = 8.2                                        # upper limit

Z1 = (math.log(b) - alpha)/float(0.1)            # Z value correponding to upper limit
Z2 = (math.log(a) - alpha)/float(0.1)            # Z value correponding to lower limit

p1 = 0.8413                                    # probability corresponding to Z1
p2 = 0.0274                                    # probability corresponding to Z2
p = round(p1 - p2,4)                           # Required probability

# Result
print " required probability: ",p

 required probability:  0.8139

Example, Page-138

In [12]:
# Variable declaration
Mean = 4.0
std_dev = 0.3

# Calculation
from scipy import *
from pylab import *

val = 33
index = (val-Mean)/std_dev

f = 0.0465     # F value for F(index) i.e. F(-1.68) = 0.9535

# Result
print "Probability: ",1-f

Probability:  0.9535

Example, page-138

In [13]:
# Variable declaration
alpha = 2               # mean of normal distribtion
beta = 0.1              # standard deviation of normal distribution

# Calculation
from scipy import *
from pylab import *

Mean = math.pow(e, 2.0 + math.pow(0.1,2.0)/2.0)                                        # mean of log normal distribution
Variance = math.pow(e, 2*2.0 + math.pow(0.1,2.0)) * ( pow(e, pow(0.1,2.0)) -1.0)        # variance of log normal distribution

Variance = round(Variance,2)
Mean = round(Mean,2)

# Resultath.
print " Mean: ",Mean
print "Variance: ",Variance

 Mean:  7.43
Variance:  0.55

Example, page-140

In [14]:
# Variable declaration
beta = 0.25

# Calculation
from scipy import *
from pylab import *
%matplotlib inline

x = array([0.894, 0.991, 0.261, 0.186, 0.311, 0.817, 2.267, 0.091, 0.139, 0.083, 
           0.235, 0.424, 0.216, 0.579, 0.429, 0.612, 0.143, 0.055, 0.752, 0.188,
           0.071, 0.159, 0.082, 1.653, 2.010, 0.258, 0.527, 1.033, 2.863, 0.365,
           0.459, 0.431, 0.092, 0.830, 1.718, 0.099, 0.162, 0.076, 0.107, 0.278,
           0.100, 0.919, 0.900, 0.093, 0.041, 0.712, 0.994, 0.149, 0.866, 0.054])
f=[]                                     # list of probability values corresponding to each x
x1=[]

for i in range(0,50):
    x1.append(x[i])

x=x1

for each in x:
    val = (1/0.25) * ( pow(e,-each/0.25) )
    f.append(val)

# Result
scatter(x,f,c='r')
bar(x,f,width=1.0,fill=False)
title("$Exponential-density-function$") 
xlabel("$decay-time(s)$")
ylabel("$density$")
Out[14]:
<matplotlib.text.Text at 0x8390860>

Example, Page-141

In [15]:
#Variable Declaration
alpha = 3
beta = 1.0/3

#Calculation
from scipy import *
from pylab import *
from scipy.integrate import *

result = quad(lambda x: 3*pow(e,-3*x), 0, 1.0/12)
result=round(result[0],3)

result1 = quad(lambda x: 3*pow(e,-3*x), 3.0/4, +inf)   # probability for 0.5 < x < +inf
result1 = round(result1[0],3)

# Result
print "probability(Part A): ",result
print "probability(Part B): ",result1

probability(Part A):  0.221
probability(Part B):  0.105

Example, Page-141

In [16]:
#Variable Declaration
alpha = 3
beta = 2

#Calculation
from scipy import *
from pylab import *
from scipy.integrate import *

Mean = 3/float(3+2)

x1 = math.gamma(3+2)
x2 = math.gamma(3)
x3 = math.gamma(2)
result = quad(lambda x: (x1/float(x2*x3))*(pow(x,float(2.0))*pow((1-x),float(1.0))), 0, 0.5)
result=round(result[0],3)

# Result
print "Part(A): ",Mean
print "Part(B): ",result # i.e. 5/16

Part(A):  0.6
Part(B):  0.313

Example, page-143

In [17]:
# Variable declaration
alpha = 0.1                                                                                          # Mean of normal distribution
beta = 0.5

# Calculation
from scipy import *
from pylab import *
from scipy.integrate import *

Mean = pow(0.1,-1.0/0.5)* math.gamma(int(1+1/0.5))                                                  # Mean life
Mean = round(Mean,0)
result = quad(lambda x: 0.1*0.5*pow(x,-0.5)*pow(e,-0.1*pow(x,0.5)) , 300, +inf)     # probability of battery life > 300 hours
result = round(result[0],3)

# Result
print "Mean lifetime: ",int(Mean),"hours"
print "probability of battery life > 300 hours: ",result

Mean lifetime:  200 hours
probability of battery life > 300 hours:  0.177

Example, Page-145

In [18]:
#Variable Declaration
x = array([[0.1,0.4,0.1],[0.2,0.2,0]])

#Calculation
from scipy import *
from pylab import *

p = x[0,2]+sum(x[1:2,1:3])

p0 = sum(x[0:2,0])  # for x1=0  
p1 = sum(x[0:2,1])  # for x1=1
p2 = sum(x[0:2,2])  # for x1=2


# Result
print "Part(A): ",p
print "Part(B): P(x1=0):",p0," P(x1=1):",p1," P(x1=2):",p2

Part(A):  0.3
Part(B): P(x1=0): 0.3  P(x1=1): 0.6  P(x1=2): 0.1

Example, Page-146

In [19]:
#Variable Declaration
x = array([[0.1,0.4,0.1],[0.2,0.2,0]])

#Calculation
from scipy import *
from pylab import *

p1 = x[1,0]/sum(x[1,0:3])  # for x1=0 & x2=1  
p2 = x[1,1]/sum(x[1,0:3])  # for x1=1 & x2=1  
p3 = x[1,2]/sum(x[1,0:3])  # for x1=2 & x2=1  

# Result
print "P(0|1):",p1," P(1|1):",p2," P(2|1):",p3
print "As P(0|1) is not equal to P(0) (0.5!=0.3) i.e. X1 & X2 are dependent" 

P(0|1): 0.5  P(1|1): 0.5  P(2|1): 0.0
As P(0|1) is not equal to P(0) (0.5!=0.3) i.e. X1 & X2 are dependent

Example, page-148

In [20]:
# Calculation
from scipy import *
from pylab import *
from scipy.integrate import *

# (A) x1 limit- (1,2)   x2 limit- (2,3)
p = dblquad(lambda x1, x2: 6*exp(-2*x1-3*x2), 2, 3, lambda x1: 1, lambda x1: 2)
p1 = round(p[0],4)

# (B) x1 limit- (0,2)   x2 limit- (2,+inf)
p = dblquad(lambda x1, x2: 6*exp(-2*x1-3*x2), 2, +inf, lambda x1: 0, lambda x1: 2)
p2 = round(p[0],4)

# Result
print "P(1<x1<2 and 2<x2<3): ",p1
print "P(0<x1<2 and 2<x2<+inf): ",p2

P(1<x1<2 and 2<x2<3):  0.0003
P(0<x1<2 and 2<x2<+inf):  0.0024

Example, Page-148

In [21]:
# Calculation
from scipy import *
from pylab import *
from scipy.integrate import *

# (A) x1 limit- (0,1)   x2 limit- (0,1)
p1 = dblquad(lambda x1, x2: 6*exp(-2*x1-3*x2), 0, 1, lambda x2: 0, lambda x2: 1)
p = round(p1[0],4)

# Result
print "P(x1<1 and x2<1): ",p

P(x1<1 and x2<1):  0.8216

Example, Page-152

In [22]:
# Variable declaration
Mean = 10       # mean: kilowatt-hours
std_dev = 3     # standard deviation: kilowatt-hours
cost = 20       # in $ per kilowatt-hours

# Calculation
from scipy import *
from pylab import *

m1 = Mean*cost
s1 = std_dev*cost
v1 = s1**2

# Result
print "Mean: ",m1,"dollar"
print "Standard Deviation: ",s1,"dollar"
print "Variance: ",v1,"dollar"

Mean:  200 dollar
Standard Deviation:  60 dollar
Variance:  3600 dollar

Example, page-154

In [23]:
# Variable declaration
Mean1 = 4                                         # Mean of X1
Variance1 = 9                                     # Variance of X1

Mean2 = -2                                        # Mean of X2
Variance2 = 5                                     # Variance of X2

# Calculation
from scipy import *
from pylab import *

# (A)  E(2*X1 + X2 - 5)
Mean = 2*Mean1 + Mean2 - 5                        # Required Mean

# (B)  Var(2*X1 + X2 - 5)
Variance = pow(2,2)*Variance1 + Variance2         # Required Variance

# Result
print "Mean of (2*X1 + X2 - 5) : ",Mean
print "Variance of (2*X1 + X2 - 5) : ",Variance

Mean of (2*X1 + X2 - 5) :  1
Variance of (2*X1 + X2 - 5) :  41

Example, page-154

In [24]:
# Variable declaration
Mean1 = 35             # Mean time in coating process( in minutes)
Variance1 = 11         # variance of time in coating process( in minutes)

Variance2 = 5          # variance of time in rinse process( in minutes)
Mean2 = 8              # Mean time in rinse process( in minutes)

# Calculation
from scipy import *
from pylab import *

Mean = Mean1 + Mean2                      # Mean of total time
Variance = Variance1 + Variance2          # Variance of total time
std_dev = int(sqrt(Variance))             # standard deviation

# Result
print "Total Mean time: ",Mean,"minutes"
print "Standard Deviation of Total time: ",std_dev,"minutes"

Total Mean time:  43 minutes
Standard Deviation of Total time:  4 minutes

Example, page-164

In [25]:
# Variable declaration
l = [67,48,76,81]                # list of observations
f=[]

# Calculation
from scipy import *
from pylab import *
%matplotlib inline

# As we know Normal scores
m1 = -0.84
m2 = -0.25
m3 = 0.25
m4 = 0.84

l = sorted(l)
f = [m1,m2,m3,m4]                # List of normal scores

# Result
scatter(f,l)
title("SimpleNormalScorePlot")
xlabel("$NormalScore$")
ylabel("$OrderedObservation$")
Out[25]:
<matplotlib.text.Text at 0xa326a58>

Example, Page-168

In [26]:
# Variable declaration
val = array([0.57, 0.74, 0.26, 0.77, 0.12])
l = []

# Calculation
from scipy import *
from pylab import *

for each in val:
    l.append(math.sqrt(-20.0*math.log(1-each,e)))
    
# Result
print "Five observations: "
for i in range(0,5):
    print round(l[i],3)," "

Five observations: 
4.108  
5.191  
2.454  
5.422  
1.599

Example, page-169

In [27]:
# Variable declaration
Mean = 50         # Mean of normal distribution
std_dev = 5       # standard deviation of normal distribution

# Calculation
from scipy import *
from pylab import *

# Computer generates 2 values 0.253 and 0.531 from uniform distribution(can be obtained by reading 3 digits at a time in TABLE-7
u1 = 0.253
u2 = 0.531

# As we know z1 = sqrt(-2*ln(u2)) * cos(2*pi*u1)  and z2 = sqrt(-2*ln(u2)) * sin(2*pi*u1)
Z1 = sqrt(-2 * (math.log(u2,e))) *( math.cos(2*pi*u1) )
Z1 = round(Z1,3)

Z2 = sqrt(-2 * (math.log(u2,e))) *( math.sin(2*pi*u1) )
Z2 = round(Z2,3)

# normal values x1 = Mean + std_dev*Z1  and x2 = Mean + std_dev*Z2
X1 = Mean + std_dev*Z1 
X2 = Mean + std_dev*Z2 

# Result
print " Normal values : ",X1," , ",X2

 Normal values :  49.895  ,  55.625