Probability in R

Summary

Recall some basic rules of probability:

• \(P(\text{nothing})=0\)
• For any even \(A\), \(0\leq P(A) \leq 1\)
• \(P(A^C)=1-P(A)\)
• If \(A\subset B\) then \(P(A)\leq P(B)\)
• For any two events \(A\) and \(B\), \(P(A\bigcup B)=P(A) + P(B)-P(A\bigcap B)\)
• For any sequence of \(n\) events \(A_1,...,A_n\): \[P\left(\bigcup_{i\in n}A_i\right)\leq \sum_{i \in n}P(A_i)\]

Recall that for conditional probability:

• \(P(A|B)=\frac{P(A\bigcap B)}{P(B)}\)
• With the multiplicative law of probability: \(P(A\bigcap B)=P(A)P(B|A)=P(B)P(B|A)\)

Bayes Rule: \[P(I|E)=\frac{P(E\bigcap I)}{P(E)}=\frac{P(I)P(E|I)}{P(I)P(E|I)+P(G)P(E|G)}\]

Probability in R: Tossing Coins and Rolling Dice

What is the probality of heads? If we’re using a fair coin, the probability should be 0.5. We can test this using R to generate data given our probability model for a fair, two-sided coin.

# Make Coin Toss Function
# -----------------------
coin_toss <- function(times=100,run=1,bias=.5){
  require(foreach)
  outcomes <- times(times) %do% {
    paste(sample(c(rep("H",len=100*bias),
                   rep("T",len=100*(1-bias))),
                 size=run,replace=T),collapse="")
  }
  return(outcomes)
}
# Flip the coin once.
coin_toss(times=1)

[1] "T"

# Flip it again.
coin_toss(times=1)

[1] "H"

# If we flip the coin 1,000 times, what do we get?
library(ggplot2)
library(dplyr)
data.frame(out=coin_toss(times=1000)) %>%
  group_by(out) %>%
  summarize(N=n()) %>%
  mutate(prop=N/sum(N)) %>%
  ggplot(aes(out,prop,label=round(prop,3))) + 
  geom_col(width=.25) +
  geom_text(vjust=1,color="white") +
  labs(x="Heads or Tails?",
       y="Proportion") +
  theme_classic() + 
  geom_hline(yintercept=0.5,linetype=3)

# Now, show how the proportion of times we get heads changes the more times we flip the coin
data.frame(prop_heads=foreach(i=1:100,.combine = "c") %do% mean(coin_toss(times=i)=="H"),
           tosses=1:100) %>% # Note the use of a for loop.
  ggplot(aes(tosses,prop_heads)) + geom_line(color="blue") +
  geom_point(color="blue") +
  geom_hline(yintercept=0.5,linetype=2) +
  theme_classic() +
  labs(x="Number of Tosses",y="Proportion Heads")

  
# Note that the probability of tails is the complement of the probability of heads
heads <- coin_toss()=="H"
mean(heads) # Probability of heads

[1] 0.46

1-mean(heads) # Probability of tails

[1] 0.54

# Conditional probability.
## What if the prability of heads was determined by the presence of a magnet (we have a special coin that is magnetic), and that if the floor is magnetic the probability of getting heads is 0.75?
## Let's assume that the magnet in the floor is activated at random. 
library(foreach)
data.frame(magnet=sample(c(0.5,0.75),size=500,replace=T)) %>%
  mutate(toss=unlist(foreach(i=magnet) %do% coin_toss(times=1,bias=i))) %>%
  group_by(magnet,toss) %>% summarize(N=n()) %>%
  mutate(prop=N/sum(N)) %>% ungroup() %>%
  mutate(magnet=car::recode(magnet,"0.50='off';0.75='on'")) %>%
  ggplot(aes(toss,prop,fill=magnet)) +
  geom_col(width=.25,position="dodge") +
  theme_classic() +
  labs(x="Heads or Tails?",y="Proportion")

NAs introduced by coercion

# So can we calculate the conditional probability of heads given the magnet is turned on?
library(broom)
tidy(glm(Heads ~ Magnet, family = binomial,
         data.frame(Magnet=sample(c(0.5,0.75),size=500,replace=T)) %>%
           mutate(Heads=unlist(foreach(i=Magnet) %do% coin_toss(times=1,bias=i))=="H",
                  Magnet=car::recode(Magnet,"0.5='Off';0.75='On'")))) %>%
  mutate(Probability_of_Heads=exp(estimate)/(1+exp(estimate)),
         Magnet=c("Off","On")) %>%
  select(Magnet,Probability_of_Heads)

NAs introduced by coercionNAs introduced by coercion

# From coins to dice...
## Create Dice Rolling Function:
dice_roll <- function(times=100,sides=6,dice=1) {
  require(foreach)
  outcome <- times(times) %do% {
    sum(sample(1:sides,size=dice,replace=T))
  }
  return(outcome)
}
ggplot(NULL,aes(dice_roll(times=10000,dice=10))) + 
  geom_histogram(fill="grey",binwidth = 1) +
  theme_classic() +
  labs(x="Sum of Values after Rolling 10 Dice",
       y="Frequency")