01-rand.Rmd

# Methods for Generating Random Variables {#rvar}

## Introduction

Most of the methods so-called *computational statistics* requires generation of random variables from specified probability distribution. In hand, we can spin wheels, roll a dice, or shuffle cards. The results are chosen randomly. However, we want the same things with computer. Here, `r`. As we know, computer cannot generate complete uniform random numbers. Instead, we generate **pseudo-random** numbers.

## Pseudo-random Numbers

```{definition, name = "Pseudo-random numbers"}
Sequence of values generated deterministically which have all the appearances of being independent $unif(0, 1)$ random variables, i.e.

$$x_1, x_2, \ldots, x_n \stackrel{iid}{\sim} unif(0, 1)$$
```


- behave *as if* following $unif(0, 1)$
- typically generated from an *initial seed*

### Linear congruential generator

<!-- Let $x_0, x_1, \ldots, x_n \in \mathbb{Z}_{+}$. -->

<!-- 1. Set $x_0$ as initial seed. -->
<!-- 2. Generate $x_n, n = 1, 2, \ldots$ recursively: -->
<!--     a. $x_n = (a x_{n - 1} + c) \mod m$ -->
<!--     b. where $a, c \in \mathbb{Z}_{+}, m: \text{modulus}$ -->
<!-- 3. Compute $u_n = \frac{x_n}{m} \in (0, 1)$ -->

Then $u_1, u_2, \ldots, u_n \sim unif(0, 1)$

\begin{algorithm}[H] \label{alg:alglcg}
  \SetAlgoLined
  \SetKwInOut{Input}{input}
  \SetKwInOut{Output}{output}
  \Input{$a, c \in \mathbb{Z}_{+}$ and modulus $m$}
  Initialize $x_0$\;
  \For{$i \leftarrow 1$ \KwTo $n$}{
    $x_i = (a x_{i - 1} + c) \mod m$\;
  }
  $u_i = \frac{x_i}{m} \in (0, 1)$\;
  \Output{$u_1, u_2, \ldots, u_n \sim unif(0, 1)$}
  \caption{Linear congruential generator}
\end{algorithm}

```{r}
lcg <- function(n, seed, a, b, m) {
  x <- rep(seed, n + 1)
  for (i in 1:n) {
    x[i + 1] <- (a * x[i] + b) %% m
  }
  x[-1] / m
}
```

```{r}
tibble(
  x = lcg(1000, 0, 1664525, 1013904223, 2^32)
) %>% 
  ggplot(aes(x = x)) +
  geom_histogram(aes(y = ..density..), bins = 30, col = gg_hcl(1))
```


### Multiplicative congruential generator

As we can expect from its name, this is congruential generator with $c = 0$.

<!-- 1. Set $x_0$ as initial seed. -->
<!-- 2. Generate $x_n, n = 1, 2, \ldots$ recursively: -->
<!--     a. $x_n = a x_{n - 1} \mod m$ -->
<!--     b. where $a \in \mathbb{Z}_{+}, m: \text{modulus}$ -->
<!-- 3. Compute $u_n = \frac{x_n}{m} \in (0, 1)$ -->

<!-- Then $u_1, u_2, \ldots, u_n \sim unif(0, 1)$ -->

\begin{algorithm}[H] \label{alg:algmcg}
  \SetAlgoLined
  \SetKwInOut{Input}{input}
  \SetKwInOut{Output}{output}
  \Input{$a, \in \mathbb{Z}_{+}$ and modulus $m$}
  Initialize $x_0$\;
  \For{$i \leftarrow 1$ \KwTo $n$}{
    $x_i = a x_{i - 1} \mod m$\;
  }
  $u_i = \frac{x_i}{m} \in (0, 1)$\;
  \Output{$u_1, u_2, \ldots, u_n \sim unif(0, 1)$}
  \caption{Multiplicative congruential generator}
\end{algorithm}

We just set `b = 0` in our `lcg()` function. The **seed must not be zero**.

```{r}
tibble(
  x = lcg(1000, 5, 1664525, 0, 2^32)
) %>% 
  ggplot(aes(x = x)) +
  geom_histogram(aes(y = ..density..), bins = 30, col = gg_hcl(1))
```

### Cycle

Generate LCG $n = 32$ with $a = 1$, $c = 1$, and $m = 16$ from the seed $x_0 = 0$.

```{r}
lcg(32, 0, 1, 1, 16)
```

Observe that we have the cycle after $m$-th number. Against this problem, we give different seed from every $(im + 1)$th random number.


## The Inverse Transform Method

```{definition, icdf, name = "Inverse of CDF"}
Since some cdf $F_X$ is not strictly increasing, we difine $F_X^{-1}(y)$ for $0 < y < 1$ by

$$F_{X}^{-1}(y) := inf \{ x : F_X(x) \ge y \}$$
```

Using this definition, we can get the following theorem.

```{theorem, probint, name = "Probability Integral Transformation"}
If $X$ is a continuous random variable with cdf $F_(x)$, then
$$U \equiv F_X(X) \sim unif(0, 1)$$
```

```{proof, name = "Probability Integral Transformation"}
Let $U \sim unif(0, 1)$. Then

\begin{equation*}
  \begin{split}
    P(F_X^{-1}(U) \le x) & = P(\inf\{t : F_X(t) = U \} \le x) \\
    & = P(U \le F_X(x)) \\
    & = F_U(F_X(x)) \\
    & = F_X(x)
  \end{split}
\end{equation*}
```

Thus, to generate $n$ random variables $\sim F_X$, we can use *uniform random numbers*.

<!-- 1. form of $F_X^{-1}(u)$ -->
<!-- 2. For each $i = 1, 2, \ldots, n$: -->
<!--     a. Generate $u_i \sim unif(0, 1)$ -->
<!--     b. $x_i = F_X^{-1}(u_i)$ -->

\begin{algorithm}[H] \label{alg:alginv1}
  \SetAlgoLined
  \SetKwInOut{Input}{input}
  \SetKwInOut{Output}{output}
  \Input{analytical form of $F_X^{-1}$}
  \For{$i \leftarrow 1$ \KwTo $n$}{
    $u_i \iid unif(0, 1)$\;
    $x_i = F_X^{-1}(u_i)$\;
  }
  \Output{$x_1, x_2, \ldots, x_n \iid F_X$}
  \caption{Inverse transformation method}
\end{algorithm}

Note that in `R`, vectorized operation would be better, i.e. generate `runif(n)` and plug it into given inverse cdf.

### Continuous case

Denote that the *probability integral transformation* holds for a continuous variable. When generating continuous random variable, applying above algorithm might work.

```{example, expon, name = "Exponential distribution"}
If $X \sim Exp(\lambda)$, then $F_X(x) = 1 - e^{-\lambda x}$. We can derive the inverse function of cdf
$$F_X^{-1}(u) = \frac{1}{\lambda}\ln(1 - u)$$
```

Note that

$$U \sim unif(0, 1) \Leftrightarrow 1 - U \sim unif(0, 1)$$

Then we just can use $U$ instead of $1 - U$.

```{r}
inv_exp <- function(n, lambda) {
  -log(runif(n)) / lambda
}
```

If we generate $x_1, \ldots, x_{500} \sim Exp(\lambda = 1)$,

```{r cdfexp, fig.cap="Inverse Transformation: Exp(1)"}
gg_curve(dexp, from = 0, to = 10) +
  geom_histogram(
    data = tibble(x = inv_exp(500, lambda = 1)),
    aes(x = x, y = ..density..),
    bins = 30,
    fill = gg_hcl(1),
    alpha = .5
  )
```


### Discrete case

<!-- 1. For each $i = 1, 2, \ldots, n$: -->
<!--     a. Generate $u_i \sim unif(0, 1)$ -->
<!--     b. Take $x_i$ s.t. $F_X(x_{i - 1}) < U \le F_X(x_i)$ -->

<!-- Collect $x_1, x_2, \ldots, x_n \sim F_X$. -->

\begin{algorithm}[H] \label{alg:alginv2}
  \SetAlgoLined
  \SetKwInOut{Input}{input}
  \SetKwInOut{Output}{output}
  \Input{analytical form of $F_X$}
  \For{$i \leftarrow 1$ \KwTo $n$}{
    $u_i \iid unif(0, 1)$\;
    Take $x_i$ s.t. $F_X(x_{i - 1}) < U \le F_X(x_i)$\;
  }
  \Output{$x_1, x_2, \ldots, x_n \iid F_X$}
  \caption{Inverse transformation method in discrete case}
\end{algorithm}

```{r, echo=FALSE}
pmf <-
  tibble(
    x = 0:4,
    p = c(.1, .2, .2, .2, .3)
  )
```

```{r, exdis, echo=FALSE}
pmf %>% 
  t() %>% 
  knitr::kable(format = "pandoc", col.names = NULL, caption = "Example of a Discrete Random Variable", longtable = TRUE)
```

```{example, dismass, name = "Discrete Random Variable"}
Consider a discrete random variable $X$ with a mass function as in Table \@ref(tab:exdis).
```

i.e.

```{r massfun, echo=FALSE, fig.cap="Probability Mass Function"}
pmf %>% 
  ggplot() +
  geom_segment(aes(x = x, y = 0, xend = x, yend = p)) +
  ylab(expression(p(x)))
```

Then we have the cdf

```{r cdfun, echo=FALSE, fig.cap="CDF of the Discrete Random Variable: Illustration for discrete case"}
pmf %>% 
  mutate(
    fx = cumsum(p),
    x_end = lead(x, default = 5),
    u = fx,
    u = ifelse(u == .5, .6, u),
    fx1 = lead(fx, default = 1),
    rand = u > fx & u <= fx1
  ) %>% 
  ggplot() +
  geom_segment(aes(x = x, y = fx, xend = x_end, yend = fx)) +
  ylab(expression(F(x))) +
  geom_segment(
    aes(x = 0, y = u, xend = x_end, yend = u, colour = rand),
    linetype = "dashed",
    arrow = arrow(length = unit(.5, "cm")),
    show.legend = FALSE
  ) +
  geom_segment(
    aes(x = x_end, y = u, xend = x_end, yend = 0, colour = rand),
    linetype = "dashed",
    arrow = arrow(length = unit(.5, "cm")),
    show.legend = FALSE
  ) +
  scale_colour_manual(
    values = c("TRUE" = gg_hcl(1), "FALSE" = "#00000000")
  )
```


Remembering the algorithm, we can implement `dplyr::case_when()` here.

```{r}
rcustom <- function(n) {
  tibble(u = runif(n)) %>% 
    mutate(
      x = case_when(
        u > 0 & u <= .1 ~ 0,
        u > .1 & u <= .3 ~ 1,
        u > .3 & u <= .5 ~ 2,
        u > .5 & u <= .7 ~ 3,
        TRUE ~ 4
      )
    ) %>% 
    select(x) %>% 
    pull()
}
```

```{r randmass, fig.cap="Generated discrete random numbers"}
tibble(x = rcustom(100)) %>% 
  count(x) %>% 
  mutate(n = n / sum(n)) %>% 
  bind_cols(px = pmf %>% select(p)) %>% # pmf table
  gather(-x, key = "key", value = "value") %>% 
  ggplot(aes(x = x, fill = key)) +
  geom_bar(aes(y = value), stat = "identity", position = "dodge", width = .2) +
  scale_fill_discrete(
    name = "Compare",
    labels = c("InvTrans", expression(p(x)))
  ) +
  ylab("prob")
```

See Figure \@ref(fig:randmass). Comparing random numbers to true pmf, the result can be said okay.

### Problems with inverse transformation

Examples \@ref(exm:expon) and \@ref(exm:dismass). We could generate these random numbers because we aware of

1. analytical $F_X$
2. $F^{-1}$

In practice, however, not all distribution have analytical $F$. Numerical computing might be possible, but it is not efficient. There are other approaches.


## The Acceptance-Rejection Method

Acceptance-rejection method does not require analytical form of cdf. What we need is our *target* density (or mass) function and *proposal* density (or mass) function. Target function is what we want to generate. Propsal function is of any random variable that is *easy to generate random numbers*. From this approach, we can generate any distribution while computation is not efficient.

|pdf or pmf|target or proposal|  
|:--------:|:--:|  
| $f$ | target|  
| $g$ | proposal - easy to generate random numbers |  

First of all, $g$ should satisfy that

$$spt f \subseteq spt g$$

Next, for some (pre-specified) $c > 0$

$$\forall x \in spt f : \frac{f(x)}{g(x)} \le c$$

<!-- For $i = 1, \ldots, n$ -->

<!-- 1. $Y \sim g(Y)$ -->
<!-- 2. $U \sim unif(0, 1) \perp\!\!\!\perp Y$ -->
<!-- 3. Accept-Reject step -->
<!--     a. Accept: $U \le \frac{f(Y)}{cg(Y)} \Rightarrow x_i = Y$ -->
<!--     b. Reject: otherwise, go to step 1 -->

<!-- Collect $x_1, x_2, \ldots, x_n \stackrel{iid}{\sim} f(x)$. -->

\begin{algorithm}[H] \label{alg:algar}
  \SetAlgoLined
  \SetKwInOut{Input}{input}
  \SetKwInOut{Output}{output}
  \Input{target $f$, proposal $g$, and $c$}
  \For{$i \leftarrow 1$ \KwTo $n$}{
    $Y \sim g(y)$\; \label{alg:argoto}
    $U \sim unif(0, 1) \ind Y$\;
    \eIf{$U \le \frac{f(Y)}{cg(Y)}$}{
      Accept $x_i = Y$\;
    }{
      go to Line \ref{alg:argoto}\;
    }
  }
  \Output{$x_1, x_2, \ldots, x_n \iid f(x)$}
  \caption{Acceptance-rejection algorithm}
\end{algorithm}

### Efficiency

```{r arprop, echo=FALSE, fig.cap="Property of AR method"}
illust_g <- function(x) {
  dbeta(x, shape1 = 3, shape2 = 2) * 2
}
#--------------------------------------
tibble(x = seq(0, 1, length.out = 101)) %>% 
  mutate(
    fx = dbeta(x, shape1 = 3, shape2 = 2),
    cgx = illust_g(x) - fx
  ) %>% 
  gather(-x, key = "density", value = "value") %>% 
  mutate(density = factor(density, levels = c("cgx", "fx"))) %>% 
  ggplot(aes(x = x)) +
  geom_area(aes(y = value, fill = density, colour = density), position = "stack", alpha = .5) +
  geom_segment(aes(x = .3, y = 0, xend = .3, yend = illust_g(.3)), linetype = "dashed") +
  geom_segment(aes(x = .3, y = 0, xend = .3, yend = dbeta(.3, shape1 = 3, shape2 = 2)), col = I("red")) +
  scale_fill_manual(
    values = c("fx" = NA, "cgx" = gg_hcl(1))
  ) +
  labs(
    x = "y",
    y = "density"
  )
```

See Figure \@ref(fig:arprop). This illustrates the motivation of A-R method. Lower one is $f(x)$ and the upper one is $cg(x)$ which covers $f$. We can see that

$$0 < \frac{f(x)}{cg(x)} \le 1$$

The algorithm takes random number from $Y \sim g$ in each recursive step $i$, which is represented as a line in the figure. At this value, the algorithm accept $Y$ as random number of $f$ if

$$U \le \frac{f(Y)}{cg(Y)}$$

Suppose that we choose a point at random on a line drawn in the figure \@ref(fig:arprop). If we get the red line, we accept. Otherwise, we reject. In other words, the *colored area is where we reject the given value*. The smaller the area is, the more efficient the algorithm will be.


```{proposition, arnote, name = "Properties of A-R Method"}
See Figure \@ref(fig:arprop).

\begin{enumerate}
  \item $\frac{f(Y)}{cg(Y)} \perp\!\!\!\perp U$
  \item $0 < \frac{f(x)}{cg(x)} \le 1$
  \item Let $N$ be the number of iterations needed to get an acceptance. Then $$N \sim Geo(p) \quad \text{where}\: p \equiv P\bigg(U \le \frac{f(Y)}{cg(Y)}\bigg)$$ and so
$$
\begin{cases}
  P(N = n) = p(1 - p)^{n - 1}I_{\{1, 2, \ldots \}}(n) \\
  E(N) = \text{average number of iterations} = \frac{1}{p}
\end{cases}
$$
  \item $X \sim Y \mid U \le \frac{f(Y)}{cg(Y)}$, i.e. $$P\bigg(Y \le y \mid U \le \frac{f(Y)}{cg(Y)}\bigg) = F_X(y)$$
\end{enumerate}
```

```{remark, areff, name = "Efficiency"}
Efficiency of the A-R method depends on $p = P\bigg(U \le \frac{f(Y)}{cg(Y)}\bigg)$. In fact,

$$E(N) = \frac{1}{p} = c$$

The algorithm becomes efficient for small $c$.
```

```{proof}
Note that

$$P\bigg( U \le \frac{f(y)}{cg(y)}, Y = y \bigg) = P\bigg(Y \le \frac{g(y)}{cg(y)} \mid Y = y \bigg)P(Y = y)$$

Since $U \sim unif(0, 1)$, $P\bigg(Y \le \frac{g(y)}{cg(y)} \mid Y = y \bigg) = \frac{f(y)}{cg(y)}$.

By construction, $P(Y = y) = g(y)$.

It follows that

\begin{equation*}
  \begin{split}
    p = P\bigg( U \le  \frac{f(y)}{cg(y)} \bigg) & = \int_{-\infty}^{\infty} P\bigg( U \le \frac{f(y)}{cg(y)}, Y = y \bigg) dy \\
    & = \int_{-\infty}^{\infty} \frac{f(y)}{cg(y)} g(y) dy \\
    & = \frac{1}{c} \int_{-\infty}^{\infty}f(y)dy \\
    & = \frac{1}{c}
  \end{split}
\end{equation*}

Hence,

$$E(N) = \frac{1}{p} = c$$

We can say that the method is efficient when the acceptance rate $p$ is large, i.e. $c$ small.
```

```{corollary, argood, name = "Efficiency of A-R Method"}
A-R method is efficient when

$g(\cdot)$ is close to $f(\cdot)$ and

have small $c$.
```

```{corollary, arc, name = "Choosing c"}
To enhance the algorithm, we might choose $c$ which satisfy

$$c = \max \bigg\{ \frac{f(x)}{g(x)} : x \in spt f \bigg\}$$
```

### Examples

```{example, arbeta, name = "Beta(a,b)"}
Let $X \sim Beta(a, b)$. Then the pdf of $X$ is given by

$$f(x) = \frac{1}{B(a, b)}x^{a - 1}(1 - x)^{b - 1}I_{(0, 1)}(x)$$
```

```{solution, arbetasol, name = "Generating Beta(a,b) with A-R method"}
Consider proposal density $g(x) = I_{(0, 1)}(x)$, i.e. $unif(0, 1)$.

To determine the optimal $c$ s.t.

$$c = \max \bigg\{ \frac{f(x)}{g(x)} : x \in (0, 1) \bigg\}$$

find the maximum of

$$\frac{f(x)}{g(x)} = \frac{1}{B(a, b)}x^{a - 1}(1 - x)^{b - 1}$$

Solve

\begin{equation*}
  \begin{split}
    \frac{d}{dx}\bigg(\frac{f(x)}{g(x)}\bigg) & = \frac{1}{B(a, b)}\Big( (a-1)x^{a-2}(1 - x)^{b - 1} - (b - 1)x^{a - 1}(1 - x)^{b - 2} \Big) \\
    & = \frac{x^{a - 2}(1 - x)^{b - 2}}{B(a, b)} \Big( (a - 1)(1 - x) - (b - 1)x \Big) \\
    & = \frac{x^{a - 2}(1 - x)^{b - 2}}{B(a, b)} \big( a - 1 - (a + b - 2)x \big) \quad = 0
  \end{split}
\end{equation*}

It follows that

$$\frac{f(x)}{g(x)} \le \frac{f(\frac{a - 1}{a + b - 2})}{g(\frac{a - 1}{a + b - 2})} = c$$

if $\frac{a - 1}{a + b - 2} \neq 0, 1$
```


```{r}
ar_beta <- function(n, a, b) {
  opt_x <- (a - 1) / (a + b - 2)
  opt_c <- dbeta(opt_x, shape1 = a, shape2 = b) / dunif(opt_x)
  X <- NULL
  N <- 0
  while (N <= n) {
    Y <- runif(n)
    U <- runif(n)
    X <- c(X, Y[U <= dbeta(Y, shape1 = a, shape2 = b) / opt_c])
    N <- length(X)
    if ( N > n ) X <- X[1:n]
  }
  X
}
```

Now we try to compare this A-R function to `R` `rbeta` function.

```{r}
gen_beta <-
  tibble(
    ar_rand = ar_beta(1000, 3, 2),
    sam = rbeta(1000, 3, 2)
  ) %>% 
  gather(key = "den", value = "value")
```


```{r betahis, fig.cap="Beta(3,2) Random numbers from each function"}
gg_curve(dbeta, from = 0, to = 1, args = list(shape1 = 3, shape2 = 2)) +
  geom_histogram(
    data = gen_beta,
    aes(x = value, y = ..density.., fill = den),
    position = "identity",
    bins = 30,
    alpha = .45
  ) +
  scale_fill_discrete(
    name = "random number",
    labels = c("AR", "rbeta")
  )
```

In the Figure \@ref(fig:betahis), the both histograms are very close to the true density curve. To see more statistically, we can draw a Q-Q plot.

```{r betaqq, fig.cap="Q-Q plot for Beta(3,2) random numbers"}
gen_beta %>% 
  ggplot(aes(sample = value)) +
  stat_qq_line(
    distribution = stats::qbeta,
    dparams = list(shape1 = 3, shape2 = 2),
    col = I("grey70"),
    size = 3.5
  ) +
  stat_qq(
    aes(colour = den),
    distribution = stats::qbeta,
    dparams = list(shape1 = 3, shape2 = 2)
  ) +
  scale_colour_discrete(
    name = "random number",
    labels = c("AR", "rbeta")
  )
```

See Figure \@ref(fig:betaqq). We have got series of numbers that are sticked to the beta distribution line.


```{example, ardiscrete, name = "A-R Method for Discrete case"}
A-R method can be also implemented to discrete case such as Example \@ref(exm:dismass).
```

```{r, exdis2, echo=FALSE}
pmf %>% 
  t() %>% 
  knitr::kable(format = "pandoc", col.names = NULL, caption = "Example of a Discrete Random Variable", longtable = TRUE)
```

```{solution, ardiscretesol, name = "Generating discrete random numbers using A-R methods"}
Consider proposal $g(x) \sim \text{Discrete unif}(0, 1, 2, 3, 4)$, i.e.

$$g(0) = g(1) = \cdots = g(4) = 0.2$$

Then we set

$$c = \max\bigg\{ \frac{p(x)}{g(x)} : x = 0, \ldots, 4 \bigg\} = \max\Big\{ 0.5, 1, 1.5 \Big\} = 1.5$$
```


## Transfomation Methods

### Continuous

```{proposition, trans1, name = "Transformation between continuous random variables"}
Relation between random variables enables generating target numbers from the others.

\begin{enumerate}
  \item $Z_1, \ldots, Z_n \iid N(0, 1) \Rightarrow \sum Z_i^2 \sim \chi^2(n)$
  \item $Y_1 \sim \chi^2(m) \ind Y_2 \sim \chi^2(n) \Rightarrow \frac{Y_1 / m}{Y_2 / n} \sim F(m, n)$
  \item $Z \sim N(0, 1) \ind Y \sim \chi^2(n) \Rightarrow \frac{Z}{\sqrt{Y / n}} \sim t(n)$
  \item $Y_1, \ldots, Y_n \iid Exp(\lambda) \Rightarrow \sum Y_i^2 Gamma(n, \lambda)$
  \item $U \sim unif(0, 1) \Rightarrow (b - a)U + a \sim unif(a, b)$
  \item $U \sim Gamma(r, \lambda) \ind V \sim Gamma(s, \lambda) \Rightarrow \frac{U}{U + V} \sim Beta(r, s)$
  \item $Z \sim N(0, 1) \Rightarrow \mu + \sigma Z \sim N(\mu, \sigma^2)$
  \item $Y \sim N(\mu, \sigma^2) \Rightarrow e^Y \sim LogNormal(\mu, \sigma^2)$
\end{enumerate}
```

```{example, transbeta, name = "Generating Beta(a, b) using rgamma"}
From Proposition \@ref(prp:trans1), we can generate $Beta(a,b)$ random numbers using $Gamma(a, 1)$ and $Gamma(b, 1)$.
```

```{r}
trans_beta <- function(n, shape1, shape2) {
  u <- rgamma(n, shape = shape1, rate = 1)
  v <- rgamma(n, shape = shape2, rate = 1)
  u / (u + v)
}
```

```{r tbetafig, echo=FALSE, fig.cap="Beta(3,2) Random numbers from each function, including transformation method"}
gen_beta2 <-
  tibble(
    ar_rand = ar_beta(1000, 3, 2),
    tran = trans_beta(1000, 3, 2),
    sam = rbeta(1000, 3, 2)
  ) %>% 
  gather(key = "den", value = "value")
#---------------------------------------
gg_curve(dbeta, from = 0, to = 1, args = list(shape1 = 3, shape2 = 2)) +
  geom_histogram(
    data = gen_beta2,
    aes(x = value, y = ..density.., fill = den),
    position = "identity",
    bins = 30,
    alpha = .3
  ) +
  scale_fill_discrete(
    name = "random number",
    labels = c("AR", "rbeta", "Transformation")
  )
```

### Box-Muller transformation

Denote that Gaussian cdf has no closed form of $F_X^{-1}$. Using polar coordiantes, we can generate Normal random numers.

```{theorem, bmnorm, name = "Box-Muller transformation"}
Let $U_1, U_2 \iid unif(0,1)$. Then

$$
\begin{cases}
  Z_1 = \sqrt{-2 \ln U_2} \cos(2\pi U_1) \\
  Z_2 = \sqrt{-2 \ln U_2} \sin(2\pi U_1)
\end{cases}
$$
```

```{proof}
Write

$$
(Z_1, Z_2)^T \sim N\bigg( \begin{bmatrix}
  0 \\
  0
\end{bmatrix}, \begin{bmatrix}
  1 & 0 \\
  0 & 1
\end{bmatrix} \bigg)
$$

Then the joint pdf is given by

$$f_{Z_1, Z_2}(x_1, x_2) = \frac{1}{2\pi}\exp\bigg(-\frac{x_1^2 + x_2^2}{2}\bigg)$$

Consider polar coordiate transformation $(R, \theta)$: $x_1 = R\cos\theta$ and $x_2 = R\sin\theta$.

Since it is also random vector,

\begin{equation*}
  \begin{split}
    f_{R, \theta}(r, \theta) & = f_{Z_1, Z_2}(x_1, x_2)\lvert J \rvert \\
    & = \frac{1}{2\pi}\exp\bigg(-\frac{x_1^2 + x_2^2}{2}\bigg)\left\lvert\begin{array}{cc}
      \frac{\partial x_1}{\partial r} & \frac{\partial x_1}{\partial \theta} \\
      \frac{\partial x_2}{\partial r} & \frac{\partial x_2}{\partial \theta}
    \end{array}\right\rvert \\
    & = \frac{1}{2\pi}\exp\bigg(-\frac{r^2}{2}\bigg)\left\lvert\begin{array}{cc}
      \frac{\partial x_1}{\partial r} & \frac{\partial x_1}{\partial \theta} \\
      \frac{\partial x_2}{\partial r} & \frac{\partial x_2}{\partial \theta}
    \end{array}\right\rvert \\
    & = \frac{r}{2\pi}\exp\bigg(-\frac{r^2}{2}\bigg)
  \end{split}
\end{equation*}

Then each marginal density function can be computed as

\begin{equation*}
  \begin{split}
    f_{\theta}(\theta) & = \int_0^\infty \frac{r}{2\pi}\exp\bigg(-\frac{r^2}{2}\bigg) dr \\
    & = \frac{1}{2\pi} I_{(0, 2\pi)}(\theta) \\
    & \stackrel{d}{=} unif(0, 2\pi)
  \end{split}
\end{equation*}

\begin{equation*}
  \begin{split}
    f_R(r) & = \int_0^\theta \frac{r}{2\pi}\exp\bigg(-\frac{r^2}{2}\bigg) d\theta \\
    & = r \exp\bigg(-\frac{r^2}{2}\bigg) I_{(0, \infty)}(r)
  \end{split}
\end{equation*}

Thus,

$$f_{R,\theta} = f_{\theta}f_R \Rightarrow R \ind \theta$$

It follows from inverse transformation theorem that

$$Z_1 = R\cos\theta = \sqrt{-2 \ln U_2} \cos(2\pi U_1)$$

and that

$$Z_2 = R\sin\theta = \sqrt{-2 \ln U_2} \sin(2\pi U_1)$$

where $U_1, U_2 \iid unif(0, 1)$
```

\begin{algorithm}[H] \label{alg:algbm}
  \SetAlgoLined
  \SetKwInOut{Input}{input}
  \SetKwInOut{Output}{output}
  \For{$i \leftarrow 1$ \KwTo $n$}{
    $U_1, U_2 \iid unif(0,1)$\;
    $z_{2i - 1} = \sqrt{-2\ln U_2}\cos(2\pi U_1)$\;
    $z_{2i} = \sqrt{-2\ln U_2}\sin(2\pi U_1)$\;
  }
  \Output{$z_1, \ldots, z_n \iid N(0,1)$}
  \caption{Box-Muller transformation}
\end{algorithm}

```{r}
bmnorm <- function(n, mean = 0, sd = 1) {
  n_bm <- ceiling(n / 2)
  tibble(
    theta = runif(n = n_bm, max = 2 * pi),
    R = sqrt(-2 * log(runif(n_bm)))
  ) %>% 
    mutate(
      x1 = R * cos(theta),
      x2 = R * sin(theta)
    ) %>% 
    gather(x1, x2, key = "key", value = "value") %>% 
    mutate(value = mean + sd * value) %>% 
    select(value) %>% 
    pull()
}
```

```{r boxnorm, fig.cap="Normal random numbers by Box-Muller transformation"}
gg_curve(dnorm, from = 0, to = 6, args = list(mean = 3, sd = 1)) +
  geom_histogram(
    data = tibble(x = bmnorm(1000, mean = 3, sd = 1)),
    aes(x = x, y = ..density..),
    bins = 30,
    fill = gg_hcl(1),
    alpha = .5
  )
```

### Discrete

```{proposition, trans2, name = "Transformation between discrete random variables"}
Relation between random variables enables generating target numbers from the others.

\begin{enumerate}
  \item $Y_1, \ldots, Y_n \iid Bernoulli(p) \Rightarrow \sum Y_i^2 \sim B(n, p)$
  \item $U \sim unif(0,1) \Rightarrow X_i = \lfloor mU \rfloor + 1$
  \item $X = \text{the number of events occurring in 1 unit of time} \sim Poisson(\lambda)$
\end{enumerate}
```

```{proposition, trans3, name = "Bernoulli process"}
Let $X_1, X_2, \ldots \iid Bernoulli(p)$.

\begin{enumerate}
  \item $N = \text{the number of trials until we see a success, i.e.} X_N = 1 \Rightarrow N \sim Geo(p)$
  \item $Y_1, \ldots, Y_r \iid Geo(p) \Rightarrow \sum\limits_{i = 1}^r Y_i = \text{the number of trials until we see r successes} \sim NegBin(r, p)$
\end{enumerate}
```

```{proposition, trans4, name = "Count process"}
Let $Y_1, Y_2, \ldots \iid Exp(\lambda)$ be interarrival times. Then

$$X = \max\{ n : \sum Y_i \le 1 \} = \text{the number of events occurring in 1 unit of time} \sim Poisson(\lambda)$$
```


## Sums and Mixtures

### Convolutions

```{definition, conv, name = "Convolution"}
Let $X_1, \ldots, X_n$ be independent and identically distributed and let $S = X_1 + \cdots X_n$. Then the distribution of $S$ is called the $n$-fold convolution of $X$ and denoted by $F_X^{*(n)}$.
```

In the last chapter, we have already seen a bunch of random variables that can be generated by summing the other.

```{example, rchi, name = "Chisquare"}
Let $Z_1, \ldots, Z_n \iid N(0, 1)$. We know from Proposition \@ref(prp:trans1) that

$$V = \sum_{i = 1}^n Z_i \sim \chi^2(n)$$
```

Building a `n` $\times$ `df` matrix can be a good strategy here. After that, `rowSums` or `colSums` ends the generation work.

```{r}
conv_chisq <- function(n, df) {
  X <-
    matrix(rnorm(n * df), nrow = n, ncol = df)^2
  rowSums(X)
}
```

```{r convchi, fig.cap="$\\chi^2$ random numbers from Normal sums"}
gg_curve(dchisq, from = 0, to = 15, args = list(df = 5)) +
  geom_histogram(
    data = tibble(x = conv_chisq(1000, df = 5)),
    aes(x = x, y = ..density..),
    bins = 30,
    fill = gg_hcl(1),
    alpha = .5
  )
```

### Mixtures

```{definition, mixprob1, name = "Discrete mixture"}
A random variable $X$ is a discrete mixture if the distribution of $X$ is a weighted sum

$$F_X(x) = \sum \theta_i F_{X_i}(x)$$

where constants $\theta_i$ are called the mixing weights or mixing probabilities.
```

```{definition, mixprob2, name = "Continuous mixture"}
A random variable $X$ is a continuous mixture if the distribution of $X$ is a weighted sum

$$F_X(x) = \int_\infty^\infty F_{X \mid Y = y} (x) f_Y(y) dy$$
```


```{example, gaussmix, name = "Mixture of several Normal distributions"}
Generate a random sample of size $1000$ from a normal location mixture with components of the mixture $N(0,1)$ and $N(3,1)$, i.e.

$$F_X = p_1 F_{X_1} + (1 - p_1) F_{X_2}$$
```

To combine samples easily, we use `foreach` library.

```{r, eval=FALSE}
library(foreach)
```

As in A-R method, Bernoullin splitting would be used.

$$
\begin{cases}
  F_{X_1} & U > p_1 \\
  F_{X_2} & \text{otherwise}
\end{cases}
$$

```{r mixturecode}
mix_norm <- function(n, p1, mean1, sd1, mean2, sd2) {
  x1 <- rnorm(n, mean = mean1, sd = sd1)
  x2 <- rnorm(n, mean = mean2, sd = sd2)
  k <- as.integer(runif(n) <= p1)
  k * x1 + (1 - k) * x2
}
```

Try various $p_1$, from 0.1 to 1. We would loop and combine by `dplyr::bind_rows()`. Reason for binding is to plot.

```{r}
mixture <-
  foreach(p1 = 0:10 / 10, .combine = bind_rows) %do% {
    tibble(
      value = mix_norm(n = 1000, p1 = p1, mean1 = 0, sd1 = 1, mean2 = 3, sd2 = 1),
      key = rep(p1, 1000)
    )
  }
```

Output is long data format. So we can easily draw a line for each group (`key`).

```{r bimodal, fig.cap="Mixture normal random number for each mixing probability"}
mixture %>% 
  ggplot(aes(x = value, colour = factor(key))) +
  stat_density(geom = "line", position = "identity") +
  scale_colour_discrete(
    name = expression(p[1]),
    labels = 0:10 / 10
  ) +
  xlab("x")
```

As $p_1$ becomes larger to $1$, the distribution becomes $N(0, 1)$. On the contrary, $p_1 = 0$ results in $N(3, 1)$.


## Multivariate Normal Random Vector

```{definition, mvn, name = "Multivariate normal random vector"}
A random vector $\mathbf{X} = (X_1, \ldots, X_p)^T$ follows multivariate normal distribution if

$$f(\mathbf{x}) = \frac{1}{(2\pi)^{\frac{p}{2}\lvert \Sigma \rvert}} \exp \bigg[ -\frac{1}{2}(\mathbf{x} \boldsymbol\mu)^T \Sigma^{-1}(\mathbf{x} \boldsymbol\mu) \bigg]$$
```

```{remark}
Let $\mathbf{Z} \sim MVN(\mathbf{0}, I)$. Then

\begin{equation}
  \Sigma^{\frac{1}{2}}\mathbf{Z} + \boldsymbol\mu \sim MVN(\boldsymbol\mu, \Sigma)
  (\#eq:stdmvn)
\end{equation}
```

From this remark, we get to generate *standard normal random vector*.

### Spectral decomposition method

Note that covariance matrix is symmetric.

```{theorem, covspec, name = "Spectral decomposition"}
Suppose that $\Sigma$ is symmetric. Then

$$\Sigma = P \Lambda P^T$$

where $(\mathbf{v}_j, \lambda_j)$ corresponding eigenvector-eigenvalue

$$
\begin{cases}
  P = \begin{bmatrix} \mathbf{v}_1 & \cdots & \mathbf{v}_p \end{bmatrix} \in \R^{p \times p} \:\text{orthogonal} \\
  \Lambda = diag(\lambda_1, \ldots, \lambda_p)
\end{cases}
$$
```

```{corollary, specsqrt}
Suppose that $\Sigma$ is symmetric. Then

$$\Sigma^{\frac{1}{2}} = P \Lambda^{\frac{1}{2}} P^T$$

where $\Lambda^{\frac{1}{2}} = diag(\sqrt{\lambda_1}, \ldots, \sqrt{\lambda_p})$
```

`eigen()` performs spectral decomposition. `$values` has eigenvalues and `$vectors` has eigenvectors. We first generate matrix that consists of standard normal random vector:

$$
\begin{bmatrix}
  Z_{11} & Z_{12} & \cdots & Z_{1p} \\
  Z_{21} & Z_{22} & \cdots & Z_{2p} \\
  \vdots & \vdots & \vdots & \vdots \\
  Z_{n1} & Z_{n2} & \cdots & Z_{np}
\end{bmatrix}
$$

Denote that each observation is row. To use Equation \@ref(eq:stdmvn), we should multiply $\Sigma^{\frac{1}{2}}$ behind this matrix, not in front of. $\boldsymbol\mu$ matrix should be also made to matrix, in form of

$$
\begin{bmatrix}
  \mu_{11} & \mu_{12} & \cdots & \mu_{1p} \\
  \mu_{11} & \mu_{22} & \cdots & \mu_{1p} \\
  \vdots & \vdots & \vdots & \vdots \\
  \mu_{11} & Z_{n2} & \cdots & \mu_{1p}
\end{bmatrix} \in \R^{n \times p}
$$

```{r}
rmvn_eigen <- function(n, mu, sig) {
  d <- length(mu)
  ev <- eigen(sig, symmetric = TRUE)
  lambda <- ev$values
  P <- ev$vectors
  sig2 <- P %*% diag(sqrt(lambda)) %*% t(P)
  Z <- matrix(rnorm(n * d), nrow = n, ncol = d)
  X <- Z %*% sig2 + matrix(mu, nrow = n, ncol = d, byrow = TRUE)
  colnames(X) <- paste0("x", 1:d)
  X %>% tbl_df()
}
```

```{r}
# mean vector -------------------------------
mu <- c(0, 1, 2)
# symmetric matrix --------------------------
sig <- matrix(numeric(9), nrow = 3, ncol = 3)
diag(sig) <- rep(1, 3)
sig[lower.tri(sig)] <- c(-.5, .5, -.5) * 2
sig <- (sig + t(sig)) / 2
```

Generate

$$\mathbf{X}_i \sim MVN\bigg((`r mu`), `r sig`\bigg)$$

```{r}
(mvn3 <- rmvn_eigen(1000, mu = mu, sig = sig))
```

```{r mvneigen, fig.cap="Multivariate normal random vector - spectral decomposition method"}
mvn3 %>% 
  GGally::ggpairs(
    lower = list(continuous = GGally::wrap(gg_scatter, size = 1))
  )
```

### Singular value decomposition

SVD can be said to be a kind of generalization of spectral decomposition. This method can be used for any matrix, i.e. non-symmetric matrix. For $\Sigma$, SVD and spectral decomposition is equivalent. However, SVD does not account for symmetric property, so this method is less efficient compared to spectral decomposition.

```{r}
rmvn_svd <- function(n, mu, sig) {
  d <- length(mu)
  S <- svd(sig)
  sig2 <- S$u %*% diag(sqrt(S$d)) %*% t(S$v)
  Z <- matrix(rnorm(n * d), nrow = n, ncol = d)
  X <- Z %*% sig2 + matrix(mu, nrow = n, ncol = d, byrow = TRUE)
  colnames(X) <- paste0("x", 1:d)
  X %>% tbl_df()
}
```

```{r mvnsvd, echo=FALSE, fig.cap="Multivariate normal random vector - svd"}
rmvn_svd(1000, mu = mu, sig = sig) %>% 
  GGally::ggpairs(
    lower = list(continuous = GGally::wrap(gg_scatter, size = 1))
  )
```

### Choleski decomposition

```{theorem, covchol, name = "Cholesky decomposition"}
Suppose that $\Sigma$ is symmetric and positive definite. Then

$$\Sigma = Q^T Q$$

where $Q$ is an upper triangular matrix.
```

```{corollary, mvnchol2}
Suppose that $\Sigma$ is symmetric and positive definite. For cholesky decomposition \@ref(thm:covchol), define

$$\Sigma^{\frac{1}{2}} = Q$$
```

`chol()` computes cholesky decomposition. In `R`, it gives upper triangular $Q$. Since some statements cholesky decomposition by $\Sigma = LL^T$ with lower triangular matrix, try not to confuse.

```{r}
rmvn_chol <- function(n, mu, sig) {
  d <- length(mu)
  sig2 <- chol(sig)
  Z <- matrix(rnorm(n * d), nrow = n, ncol = d)
  X <- Z %*% sig2 + matrix(mu, nrow = n, ncol = d, byrow = TRUE)
  colnames(X) <- paste0("x", 1:d)
  X %>% tbl_df()
}
```

```{r mvnchol, echo=FALSE, fig.cap="Multivariate normal random vector - cholesky decomposition"}
rmvn_chol(1000, mu = mu, sig = sig) %>% 
  GGally::ggpairs(
    lower = list(continuous = GGally::wrap(gg_scatter, size = 1))
  )
```


## Stochastic Processes

```{definition, stoc, name = "Stochastic process"}
A stochastic process is a collection $\{ X(t) : t \in T \}$ of random variables indexed by the set $T$. The index set $T$ could be discrete or continous.

A State space is called te set of possible values that $X(t)$ can take.
```

```{definition, dmc, name = "Discrete Time Markov Chain"}
$\{ X_n : n = 0, 1, 2, \ldots \}$ is a Discrete time markov chain on $S$ if and only if

\begin{enumerate}
  \item $S$ is at most countable
  \item Markov property $P(X_{n + 1} = j \mid X_n = i, X_{n - 1} = i_{n - 1}, \ldots, X_0 = i_0) = P(X_{n + 1} = j \mid X_n = i) = P_{ij}$
\end{enumerate}

If $P_{ij}$ is fixed, then $\{ X_n \}$ is called time homogeneous. Otherwise, it is called nonhomogeneous.
```

```{definition, rwm, name = "Random walk model"}
Let $\{ Y_n : n \in \mathbb{N} \}$ be an IID process on $S$ s.t.

$$P(Y_n = k) = p_k$$

Define

$$
S_n := \begin{cases}
  0 & n = 0 \\
  S_0 + Y_1 + \cdots + Y_n & n \in \mathbb{N}
\end{cases}
$$
```

### Gambler's ruin model

```{definition, gambler, name = "Gambler\'s ruin model"}
Let $\{ Y_n : n \in \mathbb{N} \}$ be a process on $\{ -1, 1 \}$ s.t.

$$P(Y_n = 1) = p, \quad P(Y_n = -1) = 1 - p$$

Define

$$
X_n := \begin{cases}
  a & n = 0 \\
  a + Y_1 + \cdots + Y_n & n \in \mathbb{N}
\end{cases}
$$
```

```{example, gamble, name = "Gambling with coin"}
Suppose that A and B each start with a stake of \$10, and bet \$1 on consecutive coin flips. The game ends when either one of the players has all the money. Let $S_n$ be the fortune of player A at time $n$ Then $\{ S_n, n \ge 0 \}$ is a symmetric random walk with absorbing barriers at 0 and 20. Simulate a realization of the process $\{ S_n, n \ge 0 \}$ and plot $S_n$ vs the time index from time $0$ until a barrier is reached.
```

Here we have

$$P(Y_n = 1) = P(Y_n = -1) = \frac{1}{2}$$

$$
S_n := \begin{cases}
  10 & n = 0 \\
  10 + Y_1 + \cdots + Y_n & n \in \mathbb{N}
\end{cases}
$$

```{r}
gambling <- function(begin = 10, betting = 1, prob = .5) {
  N <- begin * 2
  sa <- begin
  record <- tibble(a = begin, b = begin)
  while(all(record > 0)) {
    sa <- ifelse(runif(1) <= prob, sa + betting, sa - betting)
    record <-
      record %>% 
      bind_rows(c(a = sa, b = N - sa))
    if (sa == N) break()
  }
  record %>% 
    mutate(idx = 1:n()) %>% 
    select(idx, a, b)
}
```

```{r gambfig, fig.cap="Sample path of Gambler\'s ruin model"}
gambling(begin = 10, betting = 1, prob = .5) %>% 
  gather(-idx, key = "player", value = "fortune") %>% 
  ggplot(aes(x = idx, y = fortune, colour = player)) +
  geom_path() +
  geom_point(alpha = .5, size = 1) +
  geom_hline(yintercept = c(0, 20), col = I("grey")) +
  labs(
    x = "Betting",
    y = "Fortune"
  )
```

In Figure \@ref(fig:gambfig), we can see the result of process with $p = \frac{1}{2}$. In fact, this process with probability $0.5$ is also called *symmetric random walk*.

### Homogeneous poisson process

```{definition, coutpr, name = "Count process"}
A stochastic process $\{ N(t) : t \ge 0 \}$ where $N(t)$ is total number of events that occur by time $t$ is called counting process.

\begin{enumerate}
  \item $N(t) \ge 0$
  \item $N(t) \in \mathbb{Z}$
  \item $s \le t \Rightarrow N(s) \le N(t)$
  \item For $s < t$, $N(t) - N(s) = \:\text{the number of events that occur in}\: (s,t]$
\end{enumerate}
```

*Poisson process* is one of this counting process.

```{definition, ppdef, name = "Poisson process"}
The counting process $\{ N(t), t \ge 0 \}$ is said to be a Poisson process with rate $\lambda >0$

\begin{enumerate}
  \item $N(0) = 0$
  \item $N(t) \ind N(t + s) - N(t)$
  \item Distribution of $N(t + s) - N(t)$ is the same for all values of $t$
  \item $\lim\limits_{h \rightarrow 0} \frac{P(N(h) = 1)}{h} = \lambda$
  \item $\lim\limits_{h \rightarrow 0} \frac{P(N(h) \ge 2)}{h} = 0$
\end{enumerate}
```

```{remark}
$$\{ N(t), t \ge 0 \} \sim PP(\lambda) \Rightarrow N(t) \sim Poisson(\lambda t)$$
```

We can generate Poisson process using this relationship. However, it is slow. Thus, we find another way.

```{theorem, campbell, name = "Campbell\'s Theorem"}
Let $\{ N(t), t \ge 0 \} \sim PP(\lambda)$, let $X_t$ be the interarrival time, and let $S_{t}$ be the waiting time until $t$-th event, i.e. $S_t := \sum_{i = 1}^t X_i$. Then

$$S_1, S_2, \ldots, S_n \mid N(t) = n \stackrel{d}{=} (U_{(1)}, U_{(2)}, \ldots, U_{(n)})$$

where $U_i \sim unif(0, t)$.
```

This Campbell's theorem gives solution to the PP generation.

\begin{algorithm}[H] \label{alg:algpp}
  \SetAlgoLined
  \SetKwInOut{Input}{input}
  \SetKwInOut{Output}{output}
  \Input{end time $T$}
  Generate $N \sim Poisson(\lambda T)$\;
  For $N$, generate $U_1, \ldots, U_N \stackrel{iid}{\sim} unif(0, T)$\;
  Sort $U_1, \ldots, U_N$ in ascending order, i.e. $\{ U_{(1)}, \ldots, U_{(N)} \}$\;
  Set $S_1 = U_{(1)}, \ldots, S_N = U_{(N)}$\;
  \Output{$\begin{bmatrix} 1 & S_1 \\ 2 & S_2 \\ \cdots & \cdots \\ N & S_n \end{bmatrix}$ \\ $\Rightarrow \{ N(S_n) \} \sim PP(\lambda) \: \text{on} \: [0, T]$}
  \caption{Fast algorithm for Poisson Process}
\end{algorithm}

```{r}
rpp <- function(lambda, t0) {
  N <- rpois(1, lambda = lambda * t0)
  tibble(sn = runif(N) * t0) %>% 
    arrange(sn) %>% # arrival time
    mutate(pp = 1:n()) # N(sn) ~ PP
}
```

```{r ppfig, fig.cap="Sample path of Poisson process"}
rpp(lambda = 1, t0 = 50) %>% 
  mutate(
    true_mean = sn, # sn * lambda
    true_sd = sqrt(sn) # sn * lambda
  ) %>% 
  ggplot(aes(x = sn)) +
  geom_ribbon(
    aes(ymin = true_mean - true_sd, ymax = true_mean + true_sd), 
    fill = "grey70", 
    alpha = .5
  ) +
  geom_line(aes(y = true_mean), col = I("white"), size = 2) +
  geom_path(aes(y = pp)) +
  labs(
    x = "t",
    y = expression(N(t))
  )
```

### Nonhomogeneous poisson process

Last section, what we have seen was homogeneous PP whose distribution does not depend on $t$. On the other hand, this condition can be broken.

```{definition, npp, name = "Nonhomogeneous Poisson Process"}
The counting process $\{ N(t), t \ge 0 \}$ is said to be a Nonhomogeneous Poisson process with rate $\lambda(t) >0$ if the third condition does not hold.

\begin{enumerate}
  \item $N(0) = 0$
  \item $N(t) \ind N(t + s) - N(t)$
  \item $\lim\limits_{h \rightarrow 0} \frac{P(N(h) = 1)}{h} = \lambda$
  \item $\lim\limits_{h \rightarrow 0} \frac{P(N(h) \ge 2)}{h} = 0$
\end{enumerate}

$\lambda(t)$ is called intensity at time $t$.

$m(t) := \int_0^t \lambda(s) ds, t \ge 0$ is called mean-value function.
```

We can generate this NPP by Bernoulli splitting, which is called *thining approach*.

```{lemma, thin, name = "Thinning approach"}
Choosing $\lambda$ s.t. $\forall t \le T \lambda(t) \le \lambda$. If an event of $PP(\lambda)$ counted with

$$p(t) = \frac{\lambda(t)}{\lambda}$$

then the process follows $NPP(\lambda(t))$ on $[0, T]$
```

\begin{algorithm}[H] \label{alg:algnpp}
  \SetAlgoLined
  \SetKwInOut{Input}{input}
  \SetKwInOut{Output}{output}
  Set $t = 0, N = 0$\;
  \Repeat{$t > T$}{
    Generate $Y \sim Exp(\lambda)$\;
    Set $t \leftarrow t + Y$\;
    Generate $U \sim unif(0, 1)$\;
    \uIf{$U \le \frac{\lambda(t)}{\lambda}$}{
      Set $N \leftarrow N + 1$\;
      Set $S(N) \leftarrow t$\;
    }
  }
  \caption{Thinning algorithm}
\end{algorithm}

```{r}
npp <- function(lambda, t0, intensity) {
  t <- 0
  N <- 0
  S <- tibble(t = t, N = N)
  while (t <= t0) {
    t <- t + rexp(1, lambda)
    U <- runif(1)
    if ( U <= intensity(t) / lambda) {
      N <- N + 1
      S <-
        S %>% 
        bind_rows(tibble(t = t, N = N))
    }
  }
  S %>% slice(-1)
}
```

```{r nppfig, fig.cap="Sample path of nonhomogeneous poisson process"}
intensity <- function(x) {
  tt <- x %% 480
  rate <-
    case_when(
      tt >= 0 && tt <= 120 ~ .5,
      tt > 120 && tt <= 240 ~ 1,
      tt > 240 && tt <= 360 ~ 2,
      tt > 360 && tt <= 480 ~ 1.5
    )
  rate
}
#--------------------------
npp(2, 480, intensity = intensity) %>% 
  ggplot(aes(x = t, y = N)) +
  geom_path()
```