1.1 Algorithm

1. Derive the inverse function \(F_X^{-1}\).
   
   

2. Write a function to compute \(x = F_X^{-1}(u)\).
   
   

3. For each realization,

   1. 
      
      

   2. 
      
      

1. 
   
   
   
   
   
2. 
   
   
   
   
   

3. # write code for inverse transform example
   # f_X(x) = 3x^2, 0 <= x \<= 1

1.2 Discrete RVs

If \(X\) is a discrete random variable and \(\cdots < x_{i-1} < x_i < \cdots\) are the points of discontinuity of \(F_X(x)\), then the inverse transform is \(F_X^{-1}(u) = x_i\) where \(F_X(x_{i-1}) < u \le F_X(x_i)\). This leads to the following algorithm:

1. Generate a r.v. \(U\) from Unif\((0,1)\).

2. Select \(x_i\) where \(F_X(x_{i-1}) < U \le F_X(x_{i})\).

x <- 1:3
p <- c(0.1, 0.2, 0.7)

# write code to sample from discrete dsn
n <- 1000

2.1 Algorithm

1. Find a suitable density \(g\) and envelope \(e\).

2. Sample \(Y \sim g\).

3. Sample \(U \sim \text{Unif}(0, 1)\).

4. If \(U < f(Y)/e(Y)\), accept \(Y\).

5. Repeat from Step 2 until you have generated your desired sample size.

2.2 Envelopes

Good envelopes have the following properties:

A simple approach to finding the envelope:

Can we invert \(F(x)\) analytically?

If not, find the maximum of \(f(x)\).

# pdf function, could use dbeta() instead
f <- function(x) {
    60*x^3*(1-x)^2
}

# plot pdf
x <- seq(0, 1, length.out = 100)
ggplot() +
  geom_line(aes(x, f(x)))

envelope <- function(x) {
  ## create the envelope function
}

# Accept reject algorithm
n <- 1000 # number of samples wanted
accepted <- 0 # number of accepted samples
samples <- rep(NA, n) # store the samples here

while(accepted < n) {
  # sample y from g
  
  # sample u from uniform(0,1)
  u <- runif(1)
  
  if(u < f(y)/envelope(y)) {
    # accept
    accepted <- accepted + 1
    samples[accepted] <- y
  }
}

ggplot() +
  geom_histogram(aes(sample, y = ..density..), bins = 50, ) +
  geom_line(aes(x, f(x)), colour = "red") +
  xlab("x") + ylab("f(x)")

2.3 Why does this work?

Recall that we require \[ cg(y) \ge f(y) \qquad \forall y \in \{y: f(y) > 0\}. \] Thus,

The larger the ratio \(\frac{f(y)}{cg(y)}\), the more the random variable \(Y\) looks like a random variable distributed with pdf \(f\) and the more likely \(Y\) is to be accepted.

2.4 Additional Resources

See p.g. 69-70 of Rizzo for a proof of the validity of the method.

3.1 Algorithm

Let \(X_1, \dots, X_p\) be a set of independent random variables with pdfs \(f_{X_1}, \dots, f_{X_p}\), respectively, and let \(g(X_1, \dots, X_p)\) be some transformation we are interested in simulating from.

1. Simulate \(X_1 \sim f_{X_1}, \dots, X_p \sim f_{X_p}\).

2. Compute \(G = g(X_1, \dots, X_p)\). This is one draw from \(g(X_1, \dots, X_p)\).

3. Repeat Steps 1-2 many times to simulate from the target distribution.

library(tidyverse)

# function for squared r.v.s
squares <- function(x) x^2

sample_z <- function(n, p) {
  # store the samples
  samples <- data.frame(matrix(rnorm(n*p), nrow = n))

  samples %>% 
    mutate_all("squares") %>% # square the rvs
    rowSums() # sum over rows
}

# get samples
n <- 1000 # number of samples

# apply our function over different degrees of freedom
samples <- data.frame(chisq_2 = sample_z(n, 2),
                      chisq_5 = sample_z(n, 5),
                      chisq_10 = sample_z(n, 10),
                      chisq_100 = sample_z(n, 100))

# plot results
samples %>%
  gather(distribution, sample, everything()) %>% # make easier to plot w/ facets
  separate(distribution, into = c("dsn_name", "df")) %>% # get the df
  mutate(df = as.numeric(df)) %>% # make numeric
  mutate(pdf = dchisq(sample, df)) %>% # add density function values 
  ggplot() + # plot
  geom_histogram(aes(sample, y = ..density..)) + # samples
  geom_line(aes(sample, pdf), colour = "red") + # true pdf
  facet_wrap(~df, scales = "free")

4.1 Mixtures vs. Sums

Note that mixture distributions are not the same as the distribution of a sum of r.v.s.

\(S = \frac{1}{2}(X_1 + X_2)\)

\(Z\) such that \(f_Z(z) = 0.5f_{X_1}(z) + 0.5f_{X_2}(z)\).

n <- 1000
u <- rbinom(n, 1, 0.5)

z <- u*rnorm(n) + (1 - u)*rnorm(n, 4, 1)

ggplot() +
  geom_histogram(aes(z), bins = 50)

What about \(f_Z(z) = 0.7f_{X_1}(z) + 0.3f_{X_2}(z)\)?

4.2 Models for Count Data (refresher)

Recall that the Poisson\((\lambda)\) distribution is useful for modeling count data.

\[ f(x) = \frac{\lambda^x \exp\{-\lambda\}}{x!}, \quad x = 0, 1, 2, \dots \] Where \(X =\) number of events occuring in a fixed period of time or space.

When the mean \(\lambda\) is low, then the data consists of mostly low values (i.e. \(0, 1, 2\), etc.) and less frequently higher values.

As the mean count increases, the skewness goes away and the distribution becomes approximately normal.

With the Poisson distribution, \[E[X] = Var X = \lambda.\]

fish <- read_csv("https://stats.idre.ucla.edu/stat/data/fish.csv")

# with zeroes
ggplot(fish) + geom_histogram(aes(count), binwidth = 1)

# without zeroes
fish %>%
  filter(count > 0) %>%
  ggplot() + 
  geom_histogram(aes(count), binwidth = 1)

A zero-inflated model assumes that the zero observations have two different origins – structural and sampling zeroes.

A zero-inflated model is a mixture model because the distribution is a weighted average of the sampling model (i.e. Poisson) and a point-mass at \(0\).

For \(Y\sim ZIP(\lambda)\),

\[ Y \sim \begin{cases} 0 & \text{with probability } \pi \\ \text{Poisson}(\lambda) & \text{with probability } 1-\pi \end{cases} \] So that,

\[ Y = \qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \]

To simulate from this distribution,

n <- 1000
lambda <- 5
pi <- 0.3

u <- rbinom(n, 1, pi)
zip <- u*0 + (1-u)*rpois(n, lambda)

# zero inflated model
ggplot() + geom_histogram(aes(zip), binwidth = 1)

# Poisson(5)
ggplot() + geom_histogram(aes(rpois(n, lambda)), binwidth = 1)