Gruppeoppgaver for uke 16

# Simulates some logistic processes. 16.04.07 steffeng@math.uio.no

lambda1 <- function(i) {
  res <- 0
  if ((0 <= i) && (i <= 100)) {
    res <- i - i^2/100
  }
  res
}

mu1 <- function(i) {
  i^2/100
}

lambda2 <- function(i) {
  i
}

mu2 <- function(i) {
  i^2/50
}

# Samples the waiting-time distribution
# Note that lambda and mu are two functions,
# convenient when testing several parametrizations of the logistic process.
wait <- function(lambda, mu, i) {
  u <- runif(1)
  u <- -log(u)/(lambda(i) + mu(i))
 # print(u)
}

# First (test) version which does not discard of the sample-paths
draw.old <- function(lambda, mu, X0) {
  T <- 0
  i <- X0
  X <- 0
  k <- 1
  while ((T <= 10) && (i > 0) && (k < 10^5)) {
    w <- wait(lambda, mu, i)
    T <- T + w
    if (runif(1) > (lambda(i)/(lambda(i) + mu(i)))) {
      i <- i - 1
    }
    else {
      i <- i +1
    }
    X[k] <- i
    k <- k+1    
  }
  X <- c(X0, X)
  X
}

draw <-  function(lambda, mu, X0) {
  T <- 0
  i <- X0
  k <- 1
  X <- 0
  while ((T <= 10) && (i > 0) && (k < 10^5)) {
    w <- wait(lambda, mu, i)
    T <- T + w
    if (runif(1) > (lambda(i)/(lambda(i) + mu(i)))) {
      i <- i - 1
    }
    else {
      i <- i +1
    }
    if (T <= 10) {
      X <- i
    }
    k <- k+1    
  }
  X
}

X0 <- 50

res1 <- 0
res2 <- 0

for (i in (1:1000)) {
  res1[i] <- draw(lambda1, mu1, X0)
  res2[i] <- draw(lambda2, mu2, X0)
}

par(mfrow = c(2, 1))
plot(res1, type="l")
plot(res2, type="l")

hist(res1)
hist(res2)

summary(res1)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  34.00   46.00   49.00   49.37   52.00   66.00
summary(res2)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  31.00   44.00   49.00   49.32   54.00   74.00

sd(res1)
#[1] 4.795478
sd(res2)
#[1] 7.142692

sd(res1)^2
#[1] 22.99661
sd(res2)^2
#[1] 51.01804