Pertemuan ke-4 Conditional Probability

 In addition to regular probability, we often want to figure out how probability is affected by observing some event. For example, the NFL season is rife with possibilities. From the beginning of each season, fans start trying to figure out how likely it is that their favorite team will make the playoffs. After every game the team plays, these probabilities change based on whether they won or lost. This post won't speak to how these probabilities are updated. That's the subject for a future post on Bayesian statistics. What we will explore is the concept of conditional probability, which is the probability of seeing some event knowing that some other event has actually occurred.

Some more examples of where we might encounter such conditional probabilities:

• Inveterate bridge players like my dad would keep track of cards as they got exposed
  in the pile, for that (and the bids) provided information about the likelihoods of what hand each player had. Such card counting and conditional probabilities (what's the likelihood of each hand, given what I have seen) is one of the (frowned upon) strategies for trying to beat the casinos in blackjack and poker (see the movie 21 for a Hollywood version of real-life card counting in casinos).

• When we go to the doctor to test for a disease (say tuberculosis or HIV or even,
  more commonly, strep throat and flu), we get a yes or no answer. However, no test is perfect. A positive test still means we might not have the disease, and testing negative might mean we have it, though hopefully with very little likelihood. For us, the important thing to know is, if we tested positive (an observed event), what is the chance that we truly have the disease (an unobserved event).

• Weather forecasting is based on conditional probabilities. When the forecast says that there is a 30% chance of rain, that probability is based on all the information that the meteorologists know up until that point. It's not just a roll of the dice (though sometimes, it feels that way).

Conditional Probabilities

The flu season is rapidly approaching. Each of us have some probability of getting the flu, which can be naively computed as the number of cases of flu last year divided by the number of people potentially exposed to the flu virus in that same year. Let's call this probability P(flu).

If a person gets a flu vaccination, their chance of getting the flu should change. This would be denoted as P(flu|vaccine), and is read as "probability of getting the flu given you have been vaccinated." Because of the "been vaccinated" condition, this is a conditional probability.

So how do you compute a conditional probability? There is a basic equation that defines this:

P(A and B) = P(A | B) P(B)

which also tells us that

P(A | B) = P(A and B) / P(B)

P(A and B) is often called the joint probability of A and B, and P(A) and P(B) are often called the marginal probabilities of A and B, respectively.

Adapting the equations above to our flu example,

P(Flu | Vaccine) = P(Flu and Vaccine) / P(Vaccine)

The numerator is the probability that a person gets the vaccine and the flu; the denominator is the probability that a person gets the vaccine.

Let's look at a table of hypothetical frequencies for a population:

Condition	Vaccine	No Vaccine
Flu	A	C
No Flu	B	D

So we have:

P(Flu|Vaccine) = P(Flu and Vaccine)/P(Vaccine)

Plugging in the conditions (A, B, C, & D) from our table above:

(A/(A+B+C+D))/((A+B)/(A+B+C+D)) = A/(A+B)

Next, we will swap out the the different conditions (A B C D) with numbers so that we can calculate an answer!

Condition	Vaccine	No Vaccine
Flu	10	1000
No Flu	3000	10000

Plugging in the numbers in our new table:

P(Flu|Vaccine) = P(Flu & Vaccine)/P(Vaccine)
= 10/(3010)
= 0.0033

So this probability is the chance of getting the flu only among those who were vaccinated. We have normalized the probability of an event (getting the flu) to the conditioning event (getting vaccinated) rather than to the entire sample space.

Challenge Question: According to the table above, what is the probability of getting the flu if you weren't vaccinated P(Flu | No Vaccine)? What is the probability of getting the flu P(flu) in general?

There is another way of looking at conditional probability. If we don't know anything about event B, P(A) is the size of the light blue circle within the entire sample space (denoted by the rectangle).

If we know that the conditioning event B has happened, the probability of the event A now becomes the ratio of the light blue section to the light and dark blue section. We can then make our sample space of interest the space where event B occurs.

Statistically Independent Events

We can compare the probability of an event (A) and how it changes if we know that another event (B) has happened. How does the chance of catching flu (A) change if you're vaccinated (B)? How does a football team's chance of going to the playoffs (A) change if the quarterback is injured (B)? In both these cases, we think those chances will change.

But will the chance of the Pittsburgh Steelers beating New England Patriots (sacrilegious to some, I know) in the 4 pm game depend on the Seattle Seahawks beating the San Francisco 49ers (caveat: I'm from Seattle) during the same time? We think (and hope) not.

When knowledge of one event does not change the probability of another event happening, the two events are called statistically independent. We see a lot of things that are independent in this sense. Successive tosses of a coin are independent, or so we believe. So are successive dice rolls and slot machine plays.

Statistical independence has some mathematical consequences. It implies that

P(A|B) = P(A)

which directly implies, from the definition, that

P(A and B) = P(A)P(B)

This means that we can compute the probability of two independent events happening together by merely multiplying the individual probabilities.

Caution: You'll often find probabilities of joint events like this computed as the product of the individual events. However, this is only true if the assumption of statistical independence is valid. Often times, it is not, and so you must be careful interpreting such computations.

Let's do a little experiment in R. We'll toss two fair dice, just as we did in an earlier post, and see if the results of the two dice are independent. We first roll the dice 100,000 times, and then compute the joint distribution of the results of the rolls from the two dice.

dice <- function(no_of_rolls=1){
  x <- sample(1:6, size=no_of_rolls, replace=TRUE)
  y <- sample(1:6, size=no_of_rolls, replace=TRUE)
  return(cbind(x,y))
}

set.seed(20485)
rolls <- as.data.frame(dice(100000))

library(plyr)
freq_table <- ddply(rolls, ~x, summarize,
                    y1=sum(y==1), y2=sum(y==2), y3= sum(y==3),
                    y4 = sum(y==4), y5=sum(y==5), y6=sum(y==6))
row.names(freq_table) <- paste0('x',1:6)
prob_table <- freq_table[,-1]/100000
prob_table

##        y1      y2      y3      y4      y5      y6
## x1 0.02754 0.02886 0.02723 0.02804 0.02762 0.02820
## x2 0.02656 0.02831 0.02670 0.02803 0.02807 0.02799
## x3 0.02753 0.02876 0.02745 0.02725 0.02755 0.02751
## x4 0.02783 0.02800 0.02831 0.02715 0.02806 0.02771
## x5 0.02830 0.02836 0.02668 0.02766 0.02729 0.02849
## x6 0.02826 0.02770 0.02771 0.02747 0.02757 0.02825

Let's evaluate the probability that y=1 both with and without knowledge of x. If we don't observe x, that probability is:

> mean(rolls$y==1) 

0.16602 

If we know that x=3, then the conditional probability that y=1 given x=3 is:

> mean(rolls$y[rolls$x==3]==1)

0.165793

These results are very close.

Note: R makes it very easy to do conditional probability evaluations. In R, you can restrict yourself to those observations of y when x=3 by specifying a Boolean condition as the index of the vector, as y[x==3].

If we assumed that the results from the two dice are statistically independent, we would
have, for every pair of values i,j in 1,2,3,4,5,6:

P(x=i and y=j) = P(x=i)P(y=j)

We computed the first part earlier from prob_table.

Now for the second part.

prob_x <- table(rolls$x)/100000
prob_y <- table(rolls$y)/100000

prob_table_indep <- outer(prob_x,prob_y,'*')
row.names(prob_table_indep) <- paste0('x',1:6)
colnames(prob_table_indep) <- paste0('y',1:6)
prob_table_indep

##    
##          y1      y2      y3      y4      y5      y6
##  x1 0.02781 0.02847 0.02748 0.02774 0.02783 0.02816
##  x2 0.02750 0.02816 0.02718 0.02743 0.02753 0.02786
##  x3 0.02757 0.02823 0.02725 0.02750 0.02759 0.02792
##  x4 0.02774 0.02840 0.02741 0.02767 0.02776 0.02809
##  x5 0.02769 0.02835 0.02737 0.02762 0.02771 0.02804
##  x6 0.02772 0.02838 0.02739 0.02765 0.02774 0.02807

We see that prob_table and prob_table_indep are quite close, indicating that the rolls of the two dice are probably independent.

Statistically Testing Independence

One statistical test for testing independence of two frequency distributions (which means that for any two values of x and y, their joint probability is the product of the marginal probabilities) is the Chi-squared test. In R, this is implemented by the function chisq.test.

chisq.test(rolls$x, rolls$y)

## 
## Pearson's Chi-squared test
## 
## data:  rolls$x and rolls$y
## X-squared = 23.27, df = 25, p-value = 0.5618

We see that the p-value of this test is quite large, indicating that there is insufficient evidence to suggest that x and y are not independent.

Challenge question: If two events cannot occur together (they are mutually exclusive) can they be independent?

Bayes' Rule

Understanding how conditional probabilities change as information is acquired is part of the central dogma of the Bayesian paradigm. That paradigm is based on Bayes' theorem, which is nothing but a theorem of conditional probabilities.

P(A|B) = P(B|A)P(A)/P(B)

This provides the mathematical framework for understanding how A affects B if we know something about how B affects A.

Rearranging this formula provides a bit more insight:

P(A|B)/P(A) = P(B|A)/P(B)

In other words, how knowledge of B changes the probability of A is the same as how knowledge of A changes the probability of B, at least as a ratio. If A and B are independent, this ratio is 1.

Disease Testing

Suppose we have a test for the flu that is positive 90% of the time when tested on a flu patient (P(test + | flu) = 0.9), and is negative 95% of the time when tested on a healthy person (P(test - | no flu) = 0.95). We also know that the flu is affecting about 1% of the population (P(flu)=0.01). You go to the doctor and test positive. What is the chance that you truly have the flu?

You can answer this question directly using Bayes' theorem, but we'll tackle this a bit differently. We'll create a hypothetical population of 100,000 people, and see if we can figure this out.

flu <- sample(c('No','Yes'), size=100000, replace=TRUE, prob=c(0.99,0.01))
test <- rep(NA, 100000) # create a dummy variable first
test[flu=='No'] <- sample(c('Neg','Pos'), size=sum(flu=='No'), replace=TRUE, prob=c(0.95,0.05))
test[flu=='Yes'] <- sample(c('Neg','Pos'), size=sum(flu=='Yes'), replace=TRUE,
                          prob=c(0.1, 0.9))

In the above code we first simulate who has the flu, given on average 1% of the population gets the flu. We then find out whom among those without the flu would test positive, based on P(test - | no flu) =0.95. Recall that the when considering a conditioning event, the conditioning event is considered the sample space, and so all the laws of probability hold within that space.

In particular:

P(test + | no flu) = 1 - P(test - | no flu) = 0.05

We do a similar computation for the people with flu.

The question we are asking, what is the chance that you have the flu given that you tested positive, can then be directly answered as:

mean(flu[test=='Pos']=='Yes')

0.1567

Wow! Even though the test is pretty good, the chance that we actually have the flu even if we test positive is actually pretty small. This is because the chance of actually getting the flu is pretty small in the first place. However, if we look at how much our chance of having the flu changed with a positive test, it is quite large:

mean(flu[test=='Pos']=='Yes')/mean(flu=='Yes')

15.3152

That is, the knowledge that we tested positive increased our chance of truly having the flu 15-fold! A constant issue in medicine is if we should address the absolute increase in risk (1% to 15%) or the relative risk (15-fold) when deciding on best clinical practice.

If we calculate the probability using Bayes' theorem, we get a very similar result:

P(flu=Yes|test+) = P(test+|flu=Yes)P(flu=Yes)/P(test+)

=P(test+|flu=Yes)P(flu=Yes)/[P(test+|flu=Yes)P(flu=Yes)+P(test+|flu=No)P(flu=No)]

= 0.9 x 0.01 / (0.9 x 0.01 + 0.05 x 0.99)

= 0.1538

Conclusion

Conditional probabilities and Bayes' theorem have many everyday applications such as determining the risk of our investments, what the weather will be like this weekend, and what our medical test results mean. These concepts are central to understanding the consequences of our actions and how relationships between entities can affect outcomes. With recent increases in the amount and availability of data, understanding these concepts become essential for making informed, data-driven decisions. In this post, we reviewed how to formally look at conditional probabilities, what rules they follow, how to use those rules along with Bayes' theorem to figure out the conditional probabilities of events, and even how to "flip" them.

Pertemuan ke-3 Numerical Measures

Mean

The mean of an observation variable is a numerical measure of the central location of the data values. It is the sum of its data values divided by data count.

Hence, for a data sample of size n, its sample mean is defined as follows:

Similarly, for a data population of size N, the population mean is:

Problem

Find the mean eruption duration in the data set faithful.

Solution

We apply the mean function to compute the mean value of eruptions.

> duration = faithful$eruptions     # the eruption durations 
> mean(duration)                    # apply the mean function 
[1] 3.4878

Answer

The mean eruption duration is 3.4878 minutes.

Exercise

Find the mean eruption waiting periods in faithful.

Median

The median of an observation variable is the value at the middle when the data is sorted in ascending order. It is an ordinal measure of the central location of the data values.

Problem

Find the median of the eruption duration in the data set faithful.

Solution

We apply the median function to compute the median value of eruptions.

> duration = faithful$eruptions     # the eruption durations 
> median(duration)                  # apply the median function 
[1] 4

Answer

The median of the eruption duration is 4 minutes.

Exercise

Find the median of the eruption waiting periods in faithful.

Quartile

There are several quartiles of an observation variable. The first quartile, or lower quartile, is the value that cuts off the first 25% of the data when it is sorted in ascending order. The second quartile, or median, is the value that cuts off the first 50%. The third quartile, or upper quartile, is the value that cuts off the first 75%.

Problem

Find the quartiles of the eruption durations in the data set faithful.

Solution

We apply the quantile function to compute the quartiles of eruptions.

> duration = faithful$eruptions     # the eruption durations 
> quantile(duration)                # apply the quantile function 
    0%    25%    50%    75%   100% 
1.6000 2.1627 4.0000 4.4543 5.1000

Answer

The first, second and third quartiles of the eruption duration are 2.1627, 4.0000 and 4.4543 minutes respectively.

Exercise

Find the quartiles of the eruption waiting periods in faithful.

Note

There are several algorithms for the computation of quartiles. Details can be found in the R documentation via help(quantile).

Percentile

The n^th percentile of an observation variable is the value that cuts off the first n percent of the data values when it is sorted in ascending order.

Problem

Find the 32^nd, 57^th and 98^th percentiles of the eruption durations in the data set faithful.

Solution

We apply the quantile function to compute the percentiles of eruptions with the desired percentage ratios.

> duration = faithful$eruptions     # the eruption durations 
> quantile(duration, c(.32, .57, .98)) 
   32%    57%    98% 
2.3952 4.1330 4.9330

Answer

The 32^nd, 57^th and 98^th percentiles of the eruption duration are 2.3952, 4.1330 and 4.9330 minutes respectively.

Exercise

Find the 17^th, 43^rd, 67^th and 85^th percentiles of the eruption waiting periods in faithful.

Note

There are several algorithms for the computation of percentiles. Details can be found in the R documentation via help(quantile).

Range

The range of an observation variable is the difference of its largest and smallest data values. It is a measure of how far apart the entire data spreads in value.

Problem

Find the range of the eruption duration in the data set faithful.

Solution

We apply the max and min function to compute the largest and smallest values of eruptions, then take the difference.

> duration = faithful$eruptions     # the eruption durations 
> max(duration) − min(duration)     # apply the max and min functions 
[1] 3.5

Answer

The range of the eruption duration is 3.5 minutes.

Exercise

Find the range of the eruption waiting periods in faithful.

Interquartile Range

The interquartile range of an observation variable is the difference of its upper and lower quartiles. It is a measure of how far apart the middle portion of data spreads in value.

Problem

Find the interquartile range of eruption duration in the data set faithful.

Solution

We apply the IQR function to compute the interquartile range of eruptions.

> duration = faithful$eruptions     # the eruption durations 
> IQR(duration)                     # apply the IQR function 
[1] 2.2915

Answer

The interquartile range of eruption duration is 2.2915 minutes.

Exercise

Find the interquartile range of eruption waiting periods in faithful.

Box Plot

The box plot of an observation variable is a graphical representation based on its quartiles, as well as its smallest and largest values. It attempts to provide a visual shape of the data distribution.

Problem

Find the box plot of the eruption duration in the data set faithful.

Solution

We apply the boxplot function to produce the box plot of eruptions.

> duration = faithful$eruptions       # the eruption durations 
> boxplot(duration, horizontal=TRUE)  # horizontal box plot

Answer

The box plot of the eruption duration is:

Exercise

Find the box plot of the eruption waiting periods in faithful.

Variance

The variance is a numerical measure of how the data values is dispersed around the mean. In particular, the sample variance is defined as:

Similarly, the population variance is defined in terms of the population mean μ and population size N:

Problem

Find the variance of the eruption duration in the data set faithful.

Solution

We apply the var function to compute the variance of eruptions.

> duration = faithful$eruptions    # the eruption durations 
> var(duration)                    # apply the var function 
[1] 1.3027

Answer

The variance of the eruption duration is 1.3027.

Exercise

Find the variance of the eruption waiting periods in faithful.

Standard Deviation

The standard deviation of an observation variable is the square root of its variance.

Problem

Find the standard deviation of the eruption duration in the data set faithful.

Solution

We apply the sd function to compute the standard deviation of eruptions.

> duration = faithful$eruptions    # the eruption durations 
> sd(duration)                     # apply the sd function 
[1] 1.1414

Answer

The standard deviation of the eruption duration is 1.1414.

Exercise

Find the standard deviation of the eruption waiting periods in faithful.

Covariance

The covariance of two variables x and y in a data set measures how the two are linearly related. A positive covariance would indicate a positive linear relationship between the variables, and a negative covariance would indicate the opposite.

The sample covariance is defined in terms of the sample means as:

Similarly, the population covariance is defined in terms of the population mean μ_x, μ_y as:

Problem

Find the covariance of eruption duration and waiting time in the data set faithful. Observe if there is any linear relationship between the two variables.

Solution

We apply the cov function to compute the covariance of eruptions and waiting.

> duration = faithful$eruptions   # eruption durations 
> waiting = faithful$waiting      # the waiting period 
> cov(duration, waiting)          # apply the cov function 
[1] 13.978

Answer

The covariance of eruption duration and waiting time is about 14. It indicates a positive linear relationship between the two variables.

Correlation Coefficient

The correlation coefficient of two variables in a data set equals to their covariance divided by the product of their individual standard deviations. It is a normalized measurement of how the two are linearly related.

Formally, the sample correlation coefficient is defined by the following formula, where s_x and s_y are the sample standard deviations, and s_xy is the sample covariance.

Similarly, the population correlation coefficient is defined as follows, where σ_x and σ_y are the population standard deviations, and σ_xy is the population covariance.

If the correlation coefficient is close to 1, it would indicate that the variables are positively linearly related and the scatter plot falls almost along a straight line with positive slope. For -1, it indicates that the variables are negatively linearly related and the scatter plot almost falls along a straight line with negative slope. And for zero, it would indicate a weak linear relationship between the variables.

Problem

Find the correlation coefficient of eruption duration and waiting time in the data set faithful. Observe if there is any linear relationship between the variables.

Solution

We apply the cor function to compute the correlation coefficient of eruptions and waiting.

> duration = faithful$eruptions   # eruption durations 
> waiting = faithful$waiting      # the waiting period 
> cor(duration, waiting)          # apply the cor function 
[1] 0.90081

Answer

The correlation coefficient of eruption duration and waiting time is 0.90081. Since it is rather close to 1, we can conclude that the variables are positively linearly related.

Central Moment

The k^th central moment (or moment about the mean) of a data population is:

Similarly, the k^th central moment of a data sample is:

In particular, the second central moment of a population is its variance.

Problem

Find the third central moment of eruption duration in the data set faithful.

Solution

We apply the function moment from the e1071 package. As it is not in the core R library, the package has to be installed and loaded into the R workspace.

> library(e1071)                    # load e1071 
> duration = faithful$eruptions     # eruption durations 
> moment(duration, order=3, center=TRUE) 
[1] -0.6149

Answer

The third central moment of eruption duration is -0.6149.

Exercise

Find the third central moment of eruption waiting period in faithful.

Skewness

The skewness of a data population is defined by the following formula, where μ₂ and μ₃ are the second and third central moments.

Intuitively, the skewness is a measure of symmetry. As a rule, negative skewness indicates that the mean of the data values is less than the median, and the data distribution is left-skewed. Positive skewness would indicate that the mean of the data values is larger than the median, and the data distribution is right-skewed.

Problem

Find the skewness of eruption duration in the data set faithful.

Solution

We apply the function skewness from the e1071 package to compute the skewness coefficient of eruptions. As the package is not in the core R library, it has to be installed and loaded into the R workspace.

> library(e1071)                    # load e1071 
> duration = faithful$eruptions     # eruption durations 
> skewness(duration)                # apply the skewness function 
[1] -0.41355

Answer

The skewness of eruption duration is -0.41355. It indicates that the eruption duration distribution is skewed towards the left.

Exercise

Find the skewness of eruption waiting period in faithful.