STAT421_HW_14.Rmd

---
title: "STAT 421 Homework 14"
author: "Oliver Shanklin"
date: "April 30, 2019"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


## 1)

As the state i increases the birth rate increases linearly. Since the rate of death is a square function, the larger state i, the larger the death rate. When the state is above 25 the death rate is larger than the birth rate, which means the population is dying at a faster rate than it is growing. When i is under 25, the birth rate is larger than the death rate and the population is growing more than it is dying. This process with grow between 0 and 25, and grow the fastest at around state 10. The process will decay anywhere above state 25. It will decay the fastest when i tends to infinity.

## 2)


```{r}

t <- 0
i <- 5
lambda_i <- (1/5)*i + 1
mu_i <- (1/100)*i^2

pop_list <- list()
time_list <- list()

for( j in seq(1:10000)){
  
  t <- 0
  i <- 5
  
  q <- 1
  
  event_vec <- numeric()
  time_vec <- numeric()
  
  while(t<30){
    
    lambda_i <- (1/5)*i + 1
    mu_i <- (1/100)*i^2
    gamma_i <- lambda_i + mu_i
    
    arrival <- rexp(1, gamma_i)
    t <- t + arrival
    time_vec[q] <- t
    
    birth <- rbinom(1, 1, lambda_i/gamma_i)
    
    if (birth == 1){
      i <- i + 1
      event_vec[q] <- i
      q <- q + 1
    } else {
      i <- i - 1
      event_vec[q] <- i
      q <- q + 1
    }
    
  }
  
  pop_list[[j]] <- event_vec
  time_list[[j]] <- time_vec
  
  
}


y_lim = c(0, max( pop_list[[1]], pop_list[[2]],pop_list[[3]] ) )

x_lim = c(0,max( time_list[[1]] ) )

plot(time_list[[1]], pop_list[[1]], col = "red", type = "b", xlim = x_lim, ylim = y_lim)

points(time_list[[2]], pop_list[[2]], col = "blue", type = "b")

points(time_list[[3]], pop_list[[3]], col = "green", type = "b")






```

## 3)

The realizations in 2 do agree with the explanation in 1. Since the two rates cross at about 25 we should see the population "even" out at around 25 as t goes to infinity.

## 4)

```{r}

max_vec <- NULL

for (i in seq(1:length(pop_list))){
  
  max_vec[i] <- (max(pop_list[[i]]) > 40)
  
}

mean(max_vec)

```