# Unlike the other modules, we have been working so far, you have to download and install...
# To install this library in Jupyter notebook
import sys
!{sys.executable} -m pip install -q --upgrade pip
!{sys.executable} -m pip install -q statistics statsmodels scipy
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import math as m
import statistics
import scipy.stats as st
import statsmodels as sm
Section 1: (Overview of Probability for Machine Learning)
Section 2: (Hypothesis Testing)
Section 3: (Journey from Variance -> Covariance -> Correlation -> Regression)
likelihood of occurrence of an event
.The value of the probability of an event to happen can lie between 0 and 1 because the favorable number of outcomes can never cross the total number of outcomes.
To understand Probability, we normally start to predict the outcomes for the tossing of coins
, rolling of dice
, or drawing a card from a pack of playing cards
. Later we apply the same concepts in the domains of Artificial Intelligence and Machine Learning.
Independent Event is an event that does not impact the other. For example, rolling a dice twice and getting two sixes in a sequence (1/6 times 1/6)
Dependent event is an event that impact the other. For Example, drawing two cards from a pack of playing cards without replacement and getting two queens in a sequence. (4/52 times 3/52)
If we're only flipping the coin once, then there are only two possible outcomes in the sample space $\Omega$: it will either be H or T (using set notation, we could write this as $\Omega$ = {H, T}). Therefore: $$ P(H) = \frac{1}{2} = 0.5 $$ Equally: $$ P(T) = \frac{1}{2} = 0.5 $$
There is only one ace of spades
in the deck, so the probability of drawing it is: $$ P(\text{ace of spades}) = \frac{1}{52} \approx 0.019 $$
In contrast there are four aces
, so the probability of drawing an ace is: $$ P(\text{ace}) = \frac{4}{52} \approx 0.077 $$
P (A ∩ B) = P(A) . P(B)
Notee: To calculate the probability of intersection of dependent events, we use conditional probability (later).
two consecutive heads
from two tosses ($\Omega$ = {HH, HT, TH, TT}) is: $$ P(\text{HH}) = \frac{1}{4} = 0.25 $$three consecutive heads
from three tosses ($\Omega$ = {HHH, HHT, HTH, THH, HTT, THT, TTH, TTT}) is: $$ P(\text{HHH}) = \frac{1}{8} = 0.125 $$exactly two heads
from three tosses ($\Omega$ = {HHH, HHT, HTH, THH, HTT, THT, TTH, TTT}) is: $$ P(\text{event}) = \frac{3}{8} = 0.375 $$at least two heads
from three tosses ($\Omega$ = {HHH, HHT, HTH, THH, HTT, THT, TTH, TTT}) is: $$ P(\text{event}) = \frac{4}{8} = 0.5 $$
What is the probability of getting 5 consecutive heads from 100 tosses?
- I know the total sample space will be 2^100
- Can you tell me the count of favourable outcomes?
factorials
(e.g., $4! = 4 \times 3 \times 2 \times 1 = 24$), which feature prominently in combinatorics, to calculate probabilities instead of painstakingly determining all of the members of the sample space $\Omega$ and counting subsets within $\Omega$.This provide us with the numerator for event-probability equation from above: $$ P(\text{event}) = \frac{\text{# of outcomes of event}}{\text{# of outcomes in }\Omega} $$
In the case of coin-flipping (or any binary process with equally probable outcomes), the denominator can be calculated with $2^n$ (where $n$ is again the number of coin flips), so: $$ \frac{\text{# of outcomes of event}}{\text{# of outcomes in }\Omega} = \frac{3}{2^n} = \frac{3}{2^3} = \frac{3}{8} = 0.375 $$
Example: Use $n \choose k$ to calculate the probability of throwing exactly three heads in five coin tosses.
$$ {n \choose k} = {5 \choose 3} = \frac{5!}{3!(5 - 3)!} = \frac{5!}{(3!)(2!)} = \frac{5 \times 4 \times 3 \times 2 \times 1}{(3 \times 2 \times 1)(2 \times 1)} = \frac{120}{(6)(2)} = \frac{120}{12} = 10 $$$$P = \frac{10}{2^n} = \frac{10}{2^5} = \frac{10}{32} = 0.3125 $$
def prob_calc(n, k):
# num = m.factorial(n)/(m.factorial(k) * m.factorial(n-k))
num = m.comb(n,k)
den = 2**n
rv = num/den
return rv
# In ten tosses what is the probability of getting exactly two heads
print("prob_calc(10, 2): ", prob_calc(10, 2))
# In ten tosses what is the probability of getting exactly 5 head
print("prob_calc(10, 5): ", prob_calc(10, 5))
# Probabilities of getting 0,1,2,3,...10 heads in ten tosses
for i in range(0,11):
print (prob_calc(10,i), end=', ')
# Use list comprehension to iterate this function on values 0 to 5
probabilities = [prob_calc(10, h) for h in range(11)]
print(probabilities)
- This shows that in a ten coin toss experiment, there is less than 1% chance that we'll get all heads, while the probability of getting five heads is 24.6%
The law of large numbers states that the more experiments we run, the closer we will tend to get to the expected probability.
# The binomial() method is used to draw samples from a binomial distribution
# where n is the number of flips and p is the probability.
#np.random.seed(54) # for reproducibility
np.random.binomial(n=4, p=0.5)
# Let me increase the count of experiments to ten thousands
np.random.binomial(n=10000, p=0.5)
ne = np.array([1, 2, 3, 4, 5, 6, 7, 2**15, 2**16, 2**17, 2**18])
heads_count = [np.random.binomial(n=a, p=0.5) for a in ne]
proportion_heads = heads_count/ne
print(proportion_heads)
Note: Above are the probabilities. This comply to the Law of Large Numbers which says that "as we increase the number of coin tosses, the probability of getting heads gets closer and closer to 50%.
Marginal Probability: The probability of an event irrespective of the outcomes of other random variables, e.g. $P(A)$.
Joint Probability: Joint probability measures the probability of two more events occurring together and at the same time. The joint probability is also called the intersection of two or more events. We can represent this relation using a Venn diagram as well. For two events A and B, it is denoted by P (A ∩ B) or P(A, B)
- where the probability rule of product is used to find the probability of intersection of events. An important requirement of the rule of product is that the events should be independent. If one were to calculate the probability of an intersection of dependent events, then a different approach involving conditional probability would be needed.
Example 1:
Consider the following three events:
Suppose \begin{equation} P(A) = 0.15\end{equation}
Now it is natural that:
Example 2:
Where:
Note that the only difference between equation $(v)$ that defines Bayes' Theorem and equation $(i)$ that defines conditional probability, is the numerator on the right-hand side of the equation. In conditional probability, we see $P(A∩B)$ in the numerator, whereas in Bayes' theorem, we see $P(B|A)∗P(A)$ in the numerator.
Example 1: What is the probability that a card drawn from a deck of playing cards is a Queen, given that it is a card of Spades?
Since both the above equations has same denominator, so that can be ignored
The given message is a HAM
Example: Suppose you have five coins and you flip them one hundred times. What is the count of heads that you have got in each of the hundred experiments?
heads_count = np.random.binomial(n=5, p=0.5, size=100)
heads_count
Let us calculate the frequency of getting 0,1,2,3,4, and 5 heads
# With only one required argument returns the unique elements of an array passed
# With return_counts=True also returns the number of times each unique item appears in `ar`
unique_values, frequency = np.unique(ar=heads_count, return_counts=True)
unique_values, frequency
A probability distribution is a statistical function that describes all the possible values and likelihoods (probabilities) that a random variable can take within a given range.
Example: Let us draw a graph showing the probability distribution of getting heads as a continuation of above example.
prob = frequency/100
prob
plt.bar(unique_values, prob)
plt.grid(True)
plt.xlabel('Heads flips (out of 5 tosses)')
plt.ylabel('Event probability');
Above graph shows the probability distribution of getting heads, when you toss five coins a hundred times. So the probability of throwing two or three heads in five tosses is much higher than the probability of tossing zero or five heads.
A probability distribution is a statistical function that describes all the possible values and likelihoods that a random variable can take within a given range.
Discrete Probability Distribution Graph:
Discrete Probability Distribution Types:
- Binomial distribution
- Bernoulli distribution
- Multinomial distribution
- Poisson distribution
Continuous Probability Distribution Graph:
amount of tip
given by customers in a restaurant, which ranges from 0 to 10.0.Continuous Probability Distribution Types:
- Normal distribution
- Standard Normal distribution
- Log Normal distribution
- Exponential distribution
- Student’s T distribution
- Chi-square distribution
Note:
Gaussian
or Normal
Distribution is given by Friedrich Gauss, also known as the Bell Curve
.Where,
- μ = Mean Value
- σ = Standard Distribution of probability.
- If mean(μ) = 0 and standard deviation(σ) = 1, then this distribution is known to be normal distribution.
- x = Normal random variable
For normally distributed data:
Real World Examples:
Example:
import numpy as np
x = np.random.normal(loc=25, scale=3, size=10000)
print("Mean: ",np.mean(x))
print("Standard Deviation: ",np.std(x))
sns.displot(x)
plt.axvline(np.mean(x), color='orange')
plt.axvline(np.mean(x) - np.std(x), color='purple')
plt.axvline(np.mean(x) + np.std(x), color='purple')
plt.axvline(np.mean(x) - 2*np.std(x), color='green')
plt.axvline(np.mean(x) + 2*np.std(x), color='green')
plt.axvline(np.mean(x) - 3*np.std(x), color='red')
plt.axvline(np.mean(x) + 3*np.std(x), color='red')
plt.show()
Mean: 24.98988377299562 Standard Deviation: 2.9722754908420788
--------------------------------------------------------------------------- NameError Traceback (most recent call last) /var/folders/1t/g3ylw8h50cjdqmk5d6jh1qmm0000gn/T/ipykernel_53880/979852923.py in <module> 3 print("Mean: ",np.mean(x)) 4 print("Standard Deviation: ",np.std(x)) ----> 5 sns.displot(x) 6 plt.axvline(np.mean(x), color='orange') 7 plt.axvline(np.mean(x) - np.std(x), color='purple') NameError: name 'sns' is not defined
Example 1: Suppose there are 10000 students in my Data Science class. Generate the random result of the students with mean of 60 and standard deviation of 10. You have secured 85% marks. Calculate your Z-score and see how many students are above you? Also draw the graph of the normal distribution
# Let us first generate the random marks of ten thousand students, with a mean of 60 and standard deviation of 10
mu = 60
sigma = 10
np.random.seed(54)
x = np.random.normal(mu, sigma, 10000)
x
# Let us verify, whether the mean and std dev of above distribution `x` is 60 and 10 respectively
print("np.mean(x): ", np.mean(x))
print("np.std(x): ", np.std(x))
# You have secured 85 marks, let us now compute your z-score
xi = 85
z = (xi - np.mean(x))/np.std(x)
z = (xi - mu)/sigma
print("Your z-score: ", z)
Your Z-Score is 2.5 standard deviation above the mean.
# Let us calculate the number of students above you out of ten thousand and visualize this by drawing a graph.
a = len(np.where(x > 85)[0])
a
Out of 10000 students, only 59 are above you :)
sns.displot(x, color='green')
#ax.set_xlim(0, 100)
plt.axvline(mu, color='orange')
for i in [-3, -2, -1, 1, 2, 3]:
plt.axvline(mu+i*sigma, color='red')
plt.axvline(xi, color='purple')
plt.show()
Example 2: This is continuation of above example. Suppose your marks in Data Science are still 85%, but this time the mean of the overall result has increased to 90% with a standard deviation of 2. Calculate your Z-score and see how many students are above you? Also draw the graph of the normal distribution
# Let us first generate the random marks of ten thousand students, with a mean of 60 and standard deviation of 10
mu = 90
sigma = 2
np.random.seed(54)
x = np.random.normal(mu, sigma, 10000)
x
#Let us verify, whether the mean and std dev of above distribution x is 90 and 2 respectively
print("np.mean(x): ", np.mean(x))
print("np.std(x): ", np.std(x))
# You have secured 85 marks, let us now compute your z-score
xi = 85
z = (xi - np.mean(x))/np.std(x)
z = (xi - mu)/sigma
print("Your z-score: ", z)
Your Z-Score is 2.5 standard deviation below the mean.
#Let us calculate the number of students above you out of ten thousand and visualize this by drawing a graph.
a = len(np.where(x > 85)[0])
a
Out of 10000 students, 9934 are above you :(
sns.displot(x, color='gray')
#ax.set_xlim(0, 100)
plt.axvline(mu, color='orange')
for i in [-3, -2, -1, 1, 2, 3]:
plt.axvline(mu+i*sigma, color='red')
plt.axvline(xi, color='purple')
plt.show()
p-value
actually quantify the probability that a given observation would occur by chance alone.59
students attained a z-score
above 2.5 and only 66
attained a z-score
below -2.5.Compute p-value corresponding to Z-Score: https://www.statology.org/calculators/
Central Limit Theorem states
:Law of Large Numbers vs Central Limit Theorem:
Drawing samples from a normal distribution
# Draw random sample of size from a normal distribution.
x1 = np.random.normal(size=10000)
sns.displot(x1);
Trial 1: # Generate a random sample from above normal distribution x1
, of size 10 (without replacement). Draw corresponding histogram.
sample = np.random.choice(a=x1, size=10, replace=False)
sns.displot(sample, color='green')
plt.xlim(-1.5, 1.5)
plt.show()
With a smaller number of samples (10 in this case), the histogram is scattered all over and does not have a definite pattern. However, by increasing the sample size, the sampling distribution starts to resemble a normal distribution. This is the "Central Limit Theorem".
# Function is passed the population distribution from which to draw samples
# n is the number of samples to be drawn, while size is the size of each sample
# Function returns a list of means of all the samples
def myfunc(distr, n, size):
sample_means = []
for i in range(n):
sample = np.random.choice(distr, size=size, replace=False) #Generates a random sample from a given 1-D array
sample_means.append(sample.mean())
return sample_means
Trial 2: # Generate 10 random samples from above normal distribution x1
, of size 10 each (without replacement). Draw histogram of the mean of those ten samples.
#Pass the displot method the list of means returned by myfunc()
sns.displot(myfunc(distr=x1, n=10, size=10), color='green')
plt.xlim(-1.5, 1.5)
plt.show()
Trial 3: # Generate 100 random samples from above normal distribution x1
, of size 10 each (without replacement). Draw histogram of the mean of those hundred samples.
#Pass the displot method the list of means returned by myfunc()
sns.displot(myfunc(distr=x1, n=100, size=10), color='green')
plt.xlim(-1.5, 1.5)
plt.show()
Trial 4: # Generate 1000 random samples from above normal distribution x1
, of size 10 each (without replacement). Draw histogram of the mean of those thousand samples.
#Pass the displot method the list of means returned by myfunc()
sns.displot(myfunc(distr=x1, n=1000, size=10), color='green')
plt.xlim(-1.5, 1.5)
plt.show()
Drawing samples from a skewed distribution
# Draw positive or negitive skewed distribution based on first parameter
x2 = st.skewnorm.rvs(10, size=10000)
sns.displot(x2);
Trial 1: Generate 10 random samples from above normal distribution x2
, of size 10 each (without replacement). Draw histogram of the mean of those ten samples.
#Pass the displot method the list of means returned by myfunc()
sns.displot(myfunc(distr=x2, n=10, size=10), color='green')
plt.show()
Trial 2: # Generate 10,000 random samples from above normal distribution x1
, of size 10 each (without replacement). Draw histogram of the mean of those thousand samples.
#Pass the displot method the list of means returned by myfunc()
sns.displot(myfunc(distr=x2, n=10000, size=10), color='green', kde=True)
plt.show()
Drawing samples from a multimodal distribution
x3 = np.concatenate((np.random.normal(size=5000), np.random.normal(loc = 4.0, size=5000)))
sns.displot(x3);
Trial 1: Generate 10 random samples from above normal distribution x3
, of size 10 each (without replacement). Draw histogram of the mean of those ten samples.
#Pass the displot method the list of means returned by myfunc()
sns.displot(myfunc(distr=x3, n=10, size=10), color='green')
plt.show()
Trial 2: Generate 10,000 random samples from above normal distribution x3
, of size 10 each (without replacement). Draw histogram of the mean of those ten thousand samples.
#Pass the displot method the list of means returned by myfunc()
sns.displot(myfunc(distr=x3, n=10000, size=10), color='green', kde=True)
plt.show()
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import math as m
import statistics
import scipy.stats as st
import statsmodels as sm
1. Formulate the hypothesis
2. Gather random sample data from target population
3. Select appropriate Test
4. Execute/Perform the Test
5. Make a Decision
6. Check Errors
Example 1: A teacher believes that a certain teaching technique will influence the mean score of his students in semester exams. But he is not sure whether it will increase or decrease the mean score, which is currently 65%. To test this he apply the new teaching technique on his class of students. After the semester result, he applies the appropriate type of hypothesis test, and see if he can reject the Null Hypothesis.
Null Hypothesis ($H_0$): New teaching technique will have no effect on mean score 65
Alternative Hypothesis ($H_A$): New teaching technique will cause the mean score to be different than 65
$H_0: \mu = 65$ $\hspace{0.5cm}$ and $\hspace{0.5cm}$ $H_A: \mu \neq 65 $
Since, the alternative hypothesis has a not equal to sign, so the test will be a two tailed test or a non-directional test; because the teacher can test for effects in both direction.
Example 2: A biologist believes that a certain fertilizer will cause the plants to grow more in one month than they normally do (status-quo). The average growth rate is currently less than or equal to 5 centimeters.
Null Hypothesis ($H_0$): Fertilizer has no effect on the mean plant growth
Alternative Hypothesis ($H_A$): Fertilizer will cause mean plant growth to increase
$H_0: \mu$ ≤ 5 cm $\hspace{0.5cm}$ and $\hspace{0.5cm}$ $H_A: \mu > 5$ cm
Since, the alternative hypothesis has a greater than sign, so the test will be a single tailed test or a directional test or a right tailed test.
Example 3: A doctor believes that a new drug will be able to reduce the blood pressure in patients. The average blood pressure is currently greater than or equal to 150. To test this he/she gave that drug to a group of patients and test the effect by checking the blood pressure for one month and then perform the test.
Null Hypothesis ($H_0$): The mean blood pressure will remain same after using the drug
Alternative Hypothesis ($H_A$): The mean blood pressure will reduce after using drug
$H_0: \mu ≥ 150$ $\hspace{0.5cm}$ and $\hspace{0.5cm}$ $H_A: \mu < 150 $
Since, the alternative hypothesis has a less than sign, so the test will be a single tailed test or a directional test or a left tailed test.
Example:
Note:
If the sample data sets are drawn from a population, which is normally distributed we used parametric tests.
The Non-Parametric tests are also known as distribution free tests, and you can apply them on sample data sets, that are drawn from normally distributed, as well as from non-normally distributed populations.
One Sample tests allows us to compare the mean of a sample against a reference mean. While two sample tests allows us to comare the mean of two separate samples, which can be independent or paired.
Z-Test is used to compare means when the population standard deviation is known and the sample size is large (> 30).
T-Test is used to compare means when the population standard deviation is NOT known and the sample size is small (< 30).
A size $\alpha$ test for $H_0: \mu = \mu_0$ versus $H_A: \mu \neq \mu_0$
Reject $H_0$ if the test statistics is in the rejection region $$R = [t: \mid t \mid > t_{\alpha/2, n-1}]$$ and accepts if in the acceptance region $$A = [t: \mid t \mid \leq t_{\alpha/2, n-1}]$$
A size $\alpha$ test for $H_0: \mu \geq \mu_0$ versus $H_A: \mu < \mu_0$
Reject $H_0$ if the test statistics is in the rejection region $$R = [t: t < - t_{\alpha, n-1}]$$ and accepts if in the acceptance region $$A = [t: t \geq - t_{\alpha, n-1}]$$
The errors in hypothesis testing are nothing but our biases towards either null hypothesis or alternative hypothesis.
continuous variables
.Where: * $\bar{x}$ is the sample mean * $\mu_0$ is a reference mean, e.g., known population mean or "null hypothesis" mean * $s_{\bar{x}}$ is the sample standard error = $\frac{\sigma}{\sqrt n}$ * Use degree of freedom as n-1
Suppose we want to know that a certain species of Penguins in Pakistan is equal to 310 pounds. Test this hypothesis using a single sample T-Test at significance level of 95%
Step 1: Formulate the Hypothesis
Step 2: Data gathering
Since we cannot check the weight of all the penguins in Pakistan, so we took a random sample of 40 penguins and weigh them. The data is as under:
n = 40 # sample size
mu = 310. # mean of population (reference mean)
xbar = 300 # mean of the sample
s = 18.5 # standard deviation of sample
Step 3: Compute T-Score
#calculate standard error
sx = s/(n**.5)
#calculate t-score
t_value = (xbar - mu)/sx
t_value
Step 4: Compute p-value corresponding to T-Score: https://www.statology.org/t-score-p-value-calculator/
Step 5: Compare p-value with α (alpha)
p_value = 0.00149
alpha = 0.025
if(p_value < alpha):
print("The result is significant, so we Reject Null Hypothesis")
else:
print("Fail to Reject NUll Hypothesis")
Suppose we want to know whether or not the mean weight of two species of penguins (Adelie and Gentoo) in Pakistan is equal. Test this hypothesis using a two sample T-Test at significance level of 95%
Step 1: Formulate the Hypothesis
Step 2: Data gathering
Note: The two samples should have approximately same variance, If this assumption is not met, you should use Welch's Test instead
alpha = 0.025 # 95%, Because it is a two tailed test, therefore, the rejection region will be on both extremes .05/2
# Data about sample of Adelie species
x1bar = 300 # sample mean
n1 = 40 # sample size
s1 = 18.5 # standard deviation of sample
# Data about sample of Gantoo species
x2bar = 305 # sample mean
n2 = 38 # sample size
s2 = 16.7 # standard deviation of sample
Step 3: Compute T-Score
# pooled estimated standard error
sp = ((((n1-1)*(s1**2)) + ((n2-1)*(s2**2)))/(n1+n2-2))**0.5
sp
t_value = (x1bar-x2bar)/(sp * (((1/n1) + (1/n2))**0.5))
t_value
Step 4: Compute p-value corresponding to T-Score: https://www.statology.org/t-score-p-value-calculator/
Step 5: Compare p-value with α (alpha)
p_value = 0.21485
alpha = 0.025
if(p_value < alpha):
print("The result is significant, so we Reject Null Hypothesis")
else:
print("The result is not significant, so we Fail to Reject NUll Hypothesis")
Therefore, we donot have sufficient evidence to say that the mean weight of two species of penguins is different
We can use Student's paired-sample (a.k.a., dependent) t-test: $$ t = \frac{\bar{d} - \mu_0}{s_\bar{d}} $$ Where:
Here's an example:
exercise
containing heart rate (pulse), of 30 persons.Step 2: Data gathering
df = sns.load_dataset('exercise')
df
np.unique(df.diet, return_counts=True)
np.unique(df.time, return_counts=True)
np.unique(df.kind, return_counts=True)
rest = df[df.kind == 'rest']
rest
run = df[df.kind == 'running']
run
sns.boxplot(x='time', y='pulse', data=rest);
sns.boxplot(x='time', y='pulse', data=run);
Step 1: Formulate the Hypothesis
Step 2: Data gathering
# For simplicity, let's only consider one of the six experimental groups, say the walking, no-fat dieters:
walk_no = df[(df.diet == 'no fat') & (df.kind == 'walking')]
walk_no
(Note how participant 16 has a relatively low heart rate at all three timepoints, whereas participant 20 has a relatively high heart rate at all three timepoints.)
sns.boxplot(x='time', y='pulse', data=walk_no);
walk_no[walk_no.time == '1 min']
walk_no[walk_no.time == '1 min']['pulse']
min1 = walk_no[walk_no.time == '1 min']['pulse'].to_numpy()
min1
min15 = walk_no[walk_no.time == '15 min']['pulse'].to_numpy()
min15
(With paired samples, we can plot the values in a scatterplot, which wouldn't make any sense for independent samples, e.g.:)
sns.scatterplot(x=min1, y=min15)
plt.title('Heart rate of no-fat dieters (beats per minute)')
plt.xlabel('After 1 minute walking')
plt.ylabel('After 15 minutes walking');
Step 3+4: Compute T-Score and p-value using ttest_rel()
method of scipy.stats
st.ttest_rel(min1, min15)
Step 5: Compare p-value with α = 0.05 (alpha)
p_value = 0.02846
alpha = 0.025
if(p_value < alpha):
print("The result is significant, so we Reject Null Hypothesis")
else:
print("The result is not significant, so we Fail to Reject NUll Hypothesis")
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import math as m
import statistics
import scipy.stats as st
import statsmodels as sm
Variance
or standard deviation
tells us how much a (single) quantity varies w.r.t. its mean 'OR' how spread out the data is around the center of the distribution (the mean).Standard Deviation
$\frac{\sum(x_i-\bar{x})^2}{n}$
Example 1:
data1 = [2, 4, 6]
print("Population Variance: ", np.var(data1, ddof=0)) # use n degrees of freedom
print("Sample Variance: ", np.var(data1, ddof=1)) # use n-1 degrees of freedom
print("Population Std Dev: ", np.std(data1, ddof=0)) #((4+0+4)/3)**.5 # use n degrees of freedom
print("Sample Std Dev: ", np.std(data1, ddof=1)) #((4+0+4)/2)**.5 # use n-1 degrees of freedom
--------------------------------------------------------------------------- NameError Traceback (most recent call last) /var/folders/1t/g3ylw8h50cjdqmk5d6jh1qmm0000gn/T/ipykernel_50926/1434275318.py in <module> ----> 1 print("Population Std Dev: ", np.std(data1, ddof=0)) #((4+0+4)/3)**.5 # use n degrees of freedom 2 print("Sample Std Dev: ", np.std(data1, ddof=1)) #((4+0+4)/2)**.5 # use n-1 degrees of freedom NameError: name 'np' is not defined
Example 2:
data2 = [75, 69, 80, 70, 60, 63, 64, 69, 71]
print("Mean: ", np.mean(data2))
print("Std Dev: ", np.std(data2))
Low standard deviation implies that most values are close to the mean.
Example 3
data3 = [44, 95, 25, 60, 76, 81, 93, 84, 71, 33, 85, 81]
print("Mean: ", np.mean(data3))
print("Std Dev: ", np.std(data3))
Note: High standard deviation implies that the values are more broadly spread out.
Positive covariance
: Indicates that two variables tend to move in the same direction.Negative covariance
: Reveals that two variables tend to move in inverse directions.Zero covarince
: Indicates that two variables have no relationship between each other.
Example 1: Positive Covariance
weight = [61,62,73,74,82,86]
height = [157,168,170,181,191,185]
np.cov(height, weight,ddof=1) # sample variance along diagonal, and other values are covariances
--------------------------------------------------------------------------- NameError Traceback (most recent call last) /var/folders/1t/g3ylw8h50cjdqmk5d6jh1qmm0000gn/T/ipykernel_50926/2270568361.py in <module> 1 weight = [61,62,73,74,82,86] 2 height = [157,168,170,181,191,185] ----> 3 np.cov(height, weight,ddof=1) # sample variance along diagonal, and other values are covariances NameError: name 'np' is not defined
print("Variance of height: ", np.var(height, ddof=1)) # sample variance
print("Variance of weight: ", np.var(weight, ddof=1)) # sample variance
print("Covariance of weight, height: ", statistics.covariance(weight, height))
print("Covariance of height, weight: ", statistics.covariance(height, weight))
- Since the covariance is positive, so that means the two variables have a positive relationship
Example 2: Negative Covariance
drug = np.array([0, 1, 2, 3, 4, 5, 6, 7.]) # Drug dosage in ml
forgetness = np.array([1.86, 1.31, .62, .33, .09, -.67, -1.23, -1.37]) # Level of forgetfullness
np.cov(drug, forgetness)
print("Variance of Drug: ", np.var(drug, ddof=1)) # sample variance
print("Variance of Forgetfullness: ", np.var(forgetness, ddof=1)) # sample variance
print("Covariance of drug, forgetfulness: ", statistics.covariance(drug, forgetness))
print("Covariance of forgetfullness, drug: ", statistics.covariance(drug, forgetness))
- Since the covariance is negative, so that means the two variables have a negative relationship
The value of Covariance range from negative infinity to positive infinity, so it tells us the direction and not the magnitude of relationship between two random variables. On the other hand the value of correlation coefficient range from -1 to +1, so it tells us the direction as well as the magnitude of relationship between the two variables.
Another difference between the two is that Covariance is affected by the change in scale. If all the values of one variable are multiplied by a constant and all the values of another variable are multiplied, by a similar or different constant, then the covariance is changed. On the other hand correlation is not affected by the change in scale.
For electricity generation using a windmill, if the speed of the wind turbine increases, the generation output will increase accordingly. Thus, the variable speed and electricity output have a positive correlation here.
The correlation coefficient is used to measure the strength of the relationship between two quantitative variables.
There are several types of correlation coefficients. The correlation coefficient (developed by Karl Pearson) denoted with $r$ or $\rho$ is the most common and is defined by:
Example 1: Strong Positive Correlation
study = np.array([1, 2, 3, 4, 5, 6, 7.]) # study hours
gpa = np.array([1.0, 1.3, 2.5, 2.6, 3.5, 3.7, 4.0]) # gpa
print("Correlation Matrix: \n", np.corrcoef(study, gpa))
sns.scatterplot(x=study, y=gpa)
plt.title("Study Hours vs GPA of Students")
plt.xlabel("Study Hours")
plt.ylabel("GPA")
plt.show()
Example 2: Strong Negative Correlation
drug = np.array([0, 1, 2, 3, 4, 5, 6, 7.]) # Drug dosage in ml
forgetness = np.array([1.86, 1.31, .62, .33, .09, -.67, -1.23, -1.37]) # Level of forgetfullness
print("Correlation Matrix: \n", np.corrcoef(drug, forgetness))
sns.scatterplot(x=drug, y=forgetness)
plt.title("Clinical Trial")
plt.xlabel("Drug dosage (mL)")
plt.ylabel("Forgetfulness")
plt.show()
Correlation lets us know the magnitude and direction of relationship between two random variables ‘x’ and ‘y’, while Regression, on the other hand, predicts the value of the dependent variable based on the known value of the independent variable, assuming that the average mathematical relationship between two or more variables.
Regression can differentiate between dependent
and independent
variables, while correlation cannot.
Regression can measure causality
as well, while correlation just tell the direction of movement between two variables.
Regression can measure the strength and direction of both linear as well as non-linear relationship between two variables, while correlation can do this for linear relationship only.
The reg(x,y) and reg(y,x) are completely different, while the corr(x,y) and corr(y,x) are identical.
The graphical representation of a correlation is a single point, while a line visualizes a linear regression (Correlation produces a single statistics, while regression produces a statistical model).
Example 1:
study = np.array([1, 2, 3, 4, 5, 6, 7.]) # study hours
gpa = np.array([1.0, 1.3, 2.5, 2.6, 3.5, 3.7, 4.0]) # gpa
# plotting the regression line using regplot
sns.regplot(x=study, y=gpa, ci=None, order=1)
plt.title("Study Hours vs GPA of Students")
plt.xlabel("Study Hours")
plt.ylabel("GPA")
plt.show()
Example 2:
drug = np.array([0, 1, 2, 3, 4, 5, 6, 7.]) # Drug dosage in ml
forgetness = np.array([1.86, 1.31, .62, .33, .09, -.67, -1.23, -1.37]) # Level of forgetfullness
# plotting the regression line using regplot
sns.regplot(x=drug, y=forgetness, ci=None, order=1);
plt.title("Clinical Trial")
plt.xlabel("Drug dosage (mL)")
plt.ylabel("Forgetfulness")
plt.show();
A simple linear regression model is a mathematical equation that allows us to predict a response for a given predictor value.
y
changes with any given change of independent/feature/predictor variable x
.
(i) Let us try to understand the idea behind fitting the line. The linear regression algorithm first throws a random line, which slowly gets closer and closer to the points until it gets to the point that works.
(ii) Let say, algorithm throws a random line and then starts asking to the points, "What they want", one by one. It picks a random point which says come closer. So it listens to it, and move the line a bit closer to that point.
(iii) And it repeats the process for every point, and many a times until we get the line which is best-fit for the data points.
Slope (m)
: $=$ rise
/run
= vertical distance
/ horizontal distance
y-intercept (b)
: It tells how high the line touches the Y-axis (vertical distance)
Rotation
¶slope
. If we want to increase the slope rotate the line counter clock-wise
. If we want to decrease the slope rotate the line in clock-wise
direction. The point where the rotation happens is the point where the line interscets the y-axis and called Pivot
.Translation
¶y-intercept
. To move the line in upward direction, increase the y-intercept
, and to move the line in downward direction, decrease the y-intercept
.
Direction | Rotation | Translation |
---|---|---|
Upward | Counter-clockwise (Increased Slope) | Translate line Up (Increased y-intercept) |
Downward | Clockwise (Decreased Slope) | Translate line down (Decreased y-intercept) |
m = 2
, and y-intercept is c = 3
. Thus the linear equation of line will be:y = 2x + 3
Now here is a very important rule of machine learning, that states: Never make drastic step, which means move the line by very tiny amount, beacuse we are going to move this many times, may be 1000 of times. So we are going to pick a very small number to multiply the equation with, which is called Learning Rate
.
In machine learning and statistics, the learning rate is a tuning parameter in an optimization algorithm.
Let say we pick a very small number 0.01, now the second step is:
(Version 1)
¶Step 1
: Start with a Random line, and we perform several steps to fit this line better and better.
Step 2
: In machine learing, we repeat a process many times to make it better, so second step is to pick a large number for repetitions or epochs (An epoch is a term used in machine learning and indicates the number of passes of the entire training dataset the machine learning algorithm has completed). So let say we pick a number 10000.
Step 3
: Pick a small number or Learning Rate, let say 0.01.
Step 4
: Now repeat the process 10000 times (loop). Following process we are going to repeat:
Pick a Random point
Move line towards point using this method:
- If point above line, and to the right of y-axis
- Rotate counter-clockwise and translate up which means
- add 0.01 to slope
- add 0.01 to y-intercept
- If point above line, and to the left of y-axis
- Rotate clockwise and translate up which means
- subtract 0.01 to slope
- add 0.01 to y-intercept
- If point below line, and to the right of y-axis
- Rotate clockwise and translate down which means
- subtract 0.01 to slope
- subtract 0.01 to y-intercept
- If point below line, and to the left of y-axis
- Rotate counter-clockwise and translate down which means
- add 0.01 to slope
- subtract 0.01 to y-intercept
Step 5
: Enjoy the Best-fitted line ☺
5-.08
4.92
(Version 2)
¶In the above algorithm, we have four conditions based on which we perform the appropriate addition or subtraction operations on the slope and y-intercept values of the line. This algorithm aovids those four if...else conditions.
Step 1
: Start with a Random line, and we perform several steps to fit this line better and better.
Step 2
: In machine learing, we repeat a process many times to make it better, so second step is to pick a large number for repetitions or epochs. (An epoch is a term used in machine learning and indicates the number of passes of the entire training dataset the machine learning algorithm has completed). So let say we pick a number 1000.
Step 3
: Pick a small number or Learning Rate, let say 0.01.
Step 4
: Now repeat the process 1000 times (loop). Following process we are going to repeat:
Step 5
: Enjoy the Best-fitted line ☺
Note: This algorithm is the basic idea behind the trandional algorithms such as Gradient descent or square error algorithm that are mostly used for peforming linear regression in machine learning.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
X = np.array([1, 2, 3, 4, 5, 6, 7.]) # study hours
Y = np.array([1.0, 1.3, 2.5, 2.6, 3.5, 3.7, 4.0]) # gpa
# plotting the regression line using regplot
sns.scatterplot(x=X, y=Y)
plt.title("Study Hours vs GPA of Students")
plt.xlabel("Study Hours")
plt.ylabel("GPA")
plt.show()
# Initially let m = 0 and c = 0. Let L be our learning rate.
#This controls how much the value of m changes with each step. L could be a small value like 0.0001 for good accuracy.
m = 0
c = 0
L = 0.0001 # The learning Rate
epochs = 10000 # The number of iterations to perform gradient descent
n = float(len(X)) # Number of elements in X
# Performing Gradient Descent
# Calculate the partial derivative of the loss function with respect to m,
# and plug in the current values of x, y, m and c in it to obtain the derivative value D.
# finally update the current value of m and c
for i in range(epochs):
Y_pred = m*X + c # The current predicted value of Y
D_m = (-2/n) * sum(X * (Y - Y_pred)) # Derivative wrt m
D_c = (-2/n) * sum(Y - Y_pred) # Derivative wrt c
m = m - L * D_m # Update m
c = c - L * D_c # Update c
print (m, c)
# Draw the regression line
Y_pred = m*X + c
plt.scatter(X, Y)
plt.plot([min(X), max(X)], [min(Y_pred), max(Y_pred)], color='red') # regression line
plt.show()
Linear Least Squares
Method¶import pandas as pd
import numpy as np
from scipy import stats as st
import statistics
import statsmodels as sm
from matplotlib import pyplot as plt
import seaborn as sns
x = np.array([1, 2, 3, 4, 5, 6, 7.]) # study hours
y = np.array([1.0, 1.3, 2.5, 2.6, 3.5, 3.7, 4.0]) # gpa
# plotting the regression line using regplot
sns.scatterplot(x=x, y=y)
plt.title("Study Hours vs GPA of Students")
plt.xlabel("Study Hours")
plt.ylabel("GPA")
plt.show()
Regression Equation for this scenario:
study hours
and the outcome/dependent/response variable is gpa
. For this scenario, the line equation is:In above equation y
is the outcome/dependent/response variable and x
is the only feature/independent/predictor variable. m
is the slope of the line and c
is the y-intercept.
In this scenario, there are only two model parameters (m
is the slope of the line and c
is the y-intercept)
In Machine Learning it is a convention to represent the model parameters with Greek letters Baita. So we represent the model parameter y-intercept as Baita-zero and the model parameter slope as Baita-one:
So the above, equation can be re-written as regression equation as follows:
The $\epsilon$ term denotes error. For a given instance $i$, $\epsilon_i$ is a measure of the difference between the true $y_i$ and the model's estimate, $\hat{y}_i$. If the model predicts $y_i$ perfectly, then $\epsilon_i = \hat{y}_i - y_i = 0$.
Together, $\beta_0$ and $\beta_1$ are called the model coefficients. To create a model, we must "learn" the values of these coefficients. And once we've learned these coefficients, we can use the model to predict! Our objective is to find the parameters $\beta_0$ and $\beta_1$ that minimize $\epsilon$ across all the available data points.
Calculate $\beta_1$:
In the case of a model with a single predictor $x$, there is a fairly straightforward linear least squares formula we can use to estimate $\beta_1$: $$ \hat{\beta}_1 = \frac{\text{cov}(x,y)}{\sigma^2_x} $$
# To find cov(x,y) and var(x), use the covariance matrix
cov_mat = np.cov(x, y)
cov_mat
beta1 = cov_mat[0,1]/cov_mat[0,0]
beta1
Calculate $\beta_0$:
With $\hat{\beta}_1$ in hand, we can then rearrange the line equation ($y = \beta_0 + \beta_1 x$) to estimate $\beta_0$: $$ y = \beta_0 + \beta_1 x $$ $$ \beta_0 = y - \beta_1 x $$
We can use the mean of x and the mean of y for calculating $\hat{\beta_0}$ for all the data points: $$ \hat{\beta_0} = \bar{y} - \hat{\beta_1} \bar{x} $$
beta0 = np.mean(y) - beta1 * np.mean(x)
beta0
Fit the Line:
x = np.array([1, 2, 3, 4, 5, 6, 7.]) # study hours
y = np.array([1.0, 1.3, 2.5, 2.6, 3.5, 3.7, 4.0]) # gpa
sns.scatterplot(x=x, y=y)
xpointsofline = np.linspace(0, 8, 1000)
ypointsofline = beta0 + beta1 * xpointsofline
sns.lineplot(x=xpointsofline, y=ypointsofline)
plt.title("Study Hours vs GPA of Students")
plt.xlabel("Study Hours")
plt.ylabel("GPA")
plt.show()
Using the Model for Prediction:
study hours
we can now use the parameter estimates $\hat{\beta}_0$ and $\hat{\beta}_1$ to predict the GPA
of a student:x_i = 1.5
y_i = beta0 + beta1*x_i
y_i
x = np.array([1, 2, 3, 4, 5, 6, 7.]) # study hours
y = np.array([1.0, 1.3, 2.5, 2.6, 3.5, 3.7, 4.0]) # gpa
sns.scatterplot(x=x, y=y)
xpointsofline = np.linspace(0, 8, 1000)
ypointsofline = beta0 + beta1 * xpointsofline
sns.lineplot(x=xpointsofline, y=ypointsofline)
plt.scatter(x_i, y_i, marker='o',s=100, color='black')
plt.title("Study Hours vs GPA of Students")
plt.xlabel("Study Hours")
plt.ylabel("GPA")
plt.show()
x = np.array([0, 1, 2, 3, 4, 5, 6, 7.]) # Drug dosage in ml
y = np.array([1.86, 1.31, .62, .33, .09, -.67, -1.23, -1.37]) # Level of forgetfullness
sns.scatterplot(x=x, y=y)
plt.title("Clinical Trial")
plt.xlabel("Drug dosage (mL)")
plt.ylabel("Forgetfulness")
plt.show()
Regression Equation for this scenario:
Calculate $\beta_1$:
In the case of a model with a single predictor $x$, there is a fairly straightforward linear least squares formula we can use to estimate $\beta_1$: $$ \hat{\beta}_1 = \frac{\text{cov}(x,y)}{\sigma^2_x} $$
# To find cov(x,y) and var(x), use the covariance matrix
cov_mat = np.cov(x, y)
cov_mat
beta1 = cov_mat[0,1]/cov_mat[0,0]
beta1
Calculate $\beta_0$:
beta0 = np.mean(y) - beta1 * np.mean(x)
beta0
x = np.array([0, 1, 2, 3, 4, 5, 6, 7.]) # Drug dosage in ml
y = np.array([1.86, 1.31, .62, .33, .09, -.67, -1.23, -1.37]) # Level of forgetfullness
sns.scatterplot(x=x, y=y)
xpointsofline = np.linspace(0, 7, 1000)
ypointsofline = beta0 + beta1 * xpointsofline
sns.lineplot(x=xpointsofline, y=ypointsofline)
plt.title("Clinical Trial")
plt.xlabel("Drug dosage (mL)")
plt.ylabel("Forgetfulness")
plt.show()
In regression model terms, if we were provided with drug dosage
we can now use the parameter estimates $\hat{\beta}_0$ and $\hat{\beta}_1$ to predict the forgetfullness
of a patient:
$$ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i $$
x_i = 4
y_i = beta0 + beta1*x_i
y_i
x = np.array([0, 1, 2, 3, 4, 5, 6, 7.]) # Drug dosage in ml
y = np.array([1.86, 1.31, .62, .33, .09, -.67, -1.23, -1.37]) # Level of forgetfullness
sns.scatterplot(x=x, y=y)
xpointsofline = np.linspace(0, 7, 1000)
ypointsofline = beta0 + beta1 * xpointsofline
sns.lineplot(x=xpointsofline, y=ypointsofline)
plt.scatter(x_i, y_i, marker='o',s=100, color='black')
plt.title("Clinical Trial")
plt.xlabel("Drug dosage (mL)")
plt.ylabel("Forgetfulness")
plt.grid(True)
plt.show()
With data from female Adélie penguins, create a linear least squares model that predicts body mass with flipper length. Predict the mass of a female Adélie penguin that has a flipper length of 197mm.
penguins = sns.load_dataset('penguins')
penguins.head()
np.unique(penguins.species, return_counts=True)
adelie = penguins[penguins.species == 'Adelie']
adelie.head()
x = adelie[adelie.sex == 'Female']['flipper_length_mm'].to_numpy()
y = adelie[adelie.sex == 'Female']['body_mass_g'].to_numpy()/1000
sns.scatterplot(x=x, y=y)
plt.title("Female Adélie Penguins")
plt.xlabel("Flipper Length (mm)")
plt.ylabel("Body Mass (kg)")
plt.show();
--------------------------------------------------------------------------- NameError Traceback (most recent call last) /var/folders/1t/g3ylw8h50cjdqmk5d6jh1qmm0000gn/T/ipykernel_73570/1975260566.py in <module> ----> 1 sns.scatterplot(x=x, y=y) 2 plt.title("Female Adélie Penguins") 3 plt.xlabel("Flipper Length (mm)") 4 plt.ylabel("Body Mass (kg)") 5 plt.show(); NameError: name 'sns' is not defined
Calculate $\beta_1$:
In the case of a model with a single predictor $x$, there is a fairly straightforward linear least squares formula we can use to estimate $\beta_1$: $$ \hat{\beta}_1 = \frac{\text{cov}(x,y)}{\sigma^2_x} $$
cov_mat = np.cov(x, y)
cov_mat
beta1 = cov_mat[0,1]/cov_mat[0,0]
beta1
Calculate $\beta_0$:
With $\hat{\beta}_1$ in hand, we can then rearrange the line equation ($y = \beta_0 + \beta_1 x$) to estimate $\beta_0$: $$ y = \beta_0 + \beta_1 x $$ $$ \beta_0 = y - \beta_1 x $$
We can use the mean of x and the mean of y for calculating $\hat{\beta_0}$ for all the data points: $$ \hat{\beta_0} = \bar{y} - \hat{\beta_1} \bar{x} $$
beta0 = y.mean() - beta1*x.mean()
beta0
Fit the Line:
sns.scatterplot(x=x, y=y)
xline = np.linspace(170, 205, 1000)
yline = beta0 + beta1*xline
sns.lineplot(x=xline, y=yline, color='orange')
plt.title("Female Adélie Penguins")
plt.xlabel("Flipper Length (mm)")
plt.ylabel("Body Mass (kg)")
plt.show();
Using the Model for Prediction:
flipper length
we can now use the parameter estimates $\hat{\beta}_0$ and $\hat{\beta}_1$ to predict the body mass
of a penguin:flipper length
of a penguin is 175 mm. Can you predict its body mass
?x_i = 175
y_i = beta0 + beta1*x_i
y_i
sns.scatterplot(x=x, y=y)
xline = np.linspace(170, 205, 1000)
yline = beta0 + beta1*xline
sns.lineplot(x=xline, y=yline, color='orange')
plt.title("Female Adélie Penguins")
plt.xlabel("Flipper Length (mm)")
plt.ylabel("Body Mass (kg)")
plt.scatter(x_i, y_i, marker='o', s=100, color='purple');
plt.grid(True)
plt.show();