assign3.R

Nathan — Apr 7, 2014, 12:19 PM

#### SET UP CONSOLE, WORKSPACE AND OUTPUT ####
rm(list = ls())
cat("\f")

options(warn = -1)
graphics.off()
par(mfrow=c(2,2))

#### Variables ####
a = 0.5
w = 20
age = 10
x = age
n = 5
radix = 1e+05

#### Functions ####
Rx <- function(x) {
  if (x >= 0 & x <= 39) {
    return(0.999)
  } else if (x >= 40 & x <= 59) {
    return(0.998)
  } else if (x >= 60 & x <= 79) {
    return(0.995)
  } else if (x >= 80 & x <= w) {
    return(0.99)
  } else {
    return(0)
  }
}

tpx <- function(t, x=age) {
  return(((1 - ((x+t)/w))^a)/((1 - (x/w))^a))
}

tqx <- function(t, x=age) {
  return(1 - tpx(t, x))
}

kdeferredqx <- function(k,x) {
  return(tpx(k,x)-tpx(k+1,x))
}

midpoint <- function(f, a, b, n) {
  l = a + ((b - a)/n) * 0.5
  A = 0

  for (i in 0:(n - 1)) {
    A = A + f(l) * ((b - a)/n)
    l = l + ((b - a)/n)
  }

  return(A)
}

createLifeTable <- function(x) {
  lx <- numeric(0)
  dx <- numeric(0)
  qx <- numeric(0)
  px <- numeric(0)
  for (i in x:(w - 1)) {
    lx = c(lx, radix * tpx(i - x, x))
    qx = c(qx, tqx(1, i))
    dx = c(dx, (qx[(i - x) + 1] * lx[(i - x) + 1]))
    px = c(px, (1 - qx[(i - x) + 1]))
  }
  agecol = seq(x:(w-1))+(x-1)
  return(data.frame(agecol,lx, dx, qx, px))
}

# use survival function
ncurtate <- function(n) {

  total = 0
  for (k in 0:(n - 1)) {
    total = total + k * kdeferredqx(k,x)
  }

  total = total + n * tpx(n)
  return(total)
} 

# get answer from integrating the survival function
ncomplete <- function(n) {
  return(midpoint(tpx, 0, n, 1000))
}

# use survival function; solve by letting tpx = 0.5
# or use life table and find where radix/2 has been surpassed
median <- function() {
  i = 1
  # account for fact that they all may die at the end
  if (lxcol[w-x] > (radix/2)) {
    return(w - 1)
  } else {
    while (lxcol[i] > radix/2) {
      i = i + 1
    }
    return(i + x - 1)
  }
}

# calculate this based on a life table
mode <- function() {
   # mode of age at death is when the most people die calculate the mode based on some age x
   max = 0
   age = 0

   # process each entry in dx
   for (i in 1:(w - x)) {
    if (dxcol[i] > max) {
      max = dxcol[i]
      age = i
    }
  }
  # age corresponds to additional age, so add this to x and subtract 1 for the floor
  return(age + x - 1)
}

# get answer from integrating the survival function
stdevTx <- function() {
  Etx = midpoint(tpx, 0, (w - x), 1000)

  intFunction <- function(t) {
    return(t*tpx(t))
  }

  Etx_squared = 2 * midpoint(intFunction, 0, (w - x), 1000)
  return(sqrt((Etx_squared - Etx^2)))
}

stdevKx <- function() {

  Ekx <- function() {
    total = 0  

    for (k in 2:(w - x)) {
      total = total + lxcol[k]
    } 
    total = total/lxcol[1]
    return(total)
  }

  Ekx_squared <- function() {
    total = 0

    for (k in 2:(w - x)) {
      total = total + lxcol[k] * (2 * (k - 1) - 1) 
    }

    total = total/lxcol[1]
    return(total)
  }

  return(sqrt(Ekx_squared() - Ekx()^2))
}

# adjust the life.table provided by s years
adjust <- function(lifetable, s) {
  lx <- numeric(0)
  dx <- numeric(0)
  qx <- numeric(0)
  px <- numeric(0)
  agecol = lifetable[,1]

  for (i in x:(w - 1)) {
    if (i == x) {
      lx = c(lx, radix * tpx(i - x, x))
      # adjust qx column by appropriate amount
      qx = c(qx, (lifetable["qx"][((i - x) + 1), 1] * Rx(x)^s))
      dx = c(dx, (qx[(i - x) + 1] * lx[(i - x) + 1]))
      px = c(px, (1 - qx[(i - x) + 1]))

    } else if (x < i & i < (w - 2)) {
      lx = c(lx, (lx[(i - x)] - dx[(i - x)]))
      qx = c(qx, (lifetable["qx"][((i - x) + 1), 1] * Rx(i)^s))
      dx = c(dx, (qx[(i - x) + 1] * lx[(i - x) + 1]))
      px = c(px, (1 - qx[(i - x) + 1]))

      # doesnt apply the Rx adjustment to the (lw-1) row to ensure that all the
      # survivors die by the limiting age
    } else {
      lx = c(lx, (lx[(i - x)] - dx[(i - x)]))
      qx = c(qx, (lifetable["qx"][((i - x) + 1), 1] * Rx(i)^0))
      dx = c(dx, (qx[(i - x) + 1] * lx[(i - x) + 1]))
      px = c(px, (1 - qx[(i - x) + 1]))
    }
  }
  return(data.frame(agecol,lx, dx, qx, px))
}
#### OUTPUT ####
table1 = createLifeTable(0)
table2 = createLifeTable(x)
table3 = adjust(table2, 10)
table4 = adjust(table2, 25)

lxcol = table2[,2]
dxcol = table2[,3]
qxcol = table2[,4]
pxcol = table2[,5]

plot(tpx, 0, (w-x), ylab = "tpx", xlab = "T(x)", main = "Survival Function")
plot(seq(0:(w-x-1))-1, dxcol, type = "l", main = "Curve of Deaths", xlab = "T(x)", ylab = "Deaths")

cat("the n-year curtate expectation is:", ncurtate(n))

the n-year curtate expectation is: 4.161

cat("\nthe n-year complete expectation is:", ncomplete(n))

the n-year complete expectation is: 4.31

cat("\nthe median age at death is:", median())

the median age at death is: 18

cat("\nthe mode of the age at death is:", mode())

the mode of the age at death is: 19

cat("\nthe standard deviation of Tx) is:", stdevTx())

the standard deviation of Tx) is: 2.981

cat("\nthe standard deviation of K(x) is:", stdevKx())

the standard deviation of K(x) is: 2.914

lxcol = table3[,2]
dxcol = table3[,3]
qxcol = table3[,4]
pxcol = table3[,5]

cat("\nthe median age at death (s = 10) is:", median())

the median age at death (s = 10) is: 18

cat("\nthe mode of the age at death (s = 10) is:", mode())

the mode of the age at death (s = 10) is: 19

lxcol = table4[,2]
dxcol = table4[,3]
qxcol = table4[,4]
pxcol = table4[,5]

cat("\nthe median age at death (s = 25) is:", median())

the median age at death (s = 25) is: 18

cat("\nthe mode of the age at death (s = 25) is:", mode())

the mode of the age at death (s = 25) is: 19