graphics.off(); rm(list=ls()); cat("\f")
global = list(a=1.2,w=100,x=25,n=40,radix=1e+05)
Rx <- function(x, w=global$w) { # reduction factors
ifelse(x<0, 0,
ifelse(x<=39, 0.999,
ifelse(x<=59, 0.998,
ifelse(x<=79, 0.995,
ifelse(x<w, 0.990, 0)))))
}
tpx <- function(t, x=global$x, w=global$w, a=global$a) { # survival function
((1 - ((x+t)/w))/(1 - (x/w)))^a
}
tqx <- function(t, x=global$x) { # 1 - survival function
1 - tpx(t,x)
}
kdeferredqx <- function(k, x=global$x) { # prob(death in next year)
tpx(k,x) - tpx(k+1,x)
}
createLifeTable <- function(x, w=global$w, radix=global$radix) { # data.frame lifeable
data.frame(x = seq(x, w-1),
lx=radix*tpx(seq(0,w-x-1),x),
dx = radix*c(tqx(1,x), tpx(seq(0,w-x-2),x)*kdeferredqx(1,seq(x,w-1)))[1:(w-x)],
qx=tqx(1, seq(x,w-1)),
px=tpx(1, seq(x,w-1)))
}
ncurtate <- function(n, x=global$x) { # curtate expected lifetime
sum(0:(n-1)*kdeferredqx(0:(n-1), x)) + n*tpx(n,x)
}
ncomplete <- function(n) { # complete expected lifetime
integrate(tpx, 0, n)$value
}
median <- function(x, w=global$w, radix=global$radix, s=0) { # median from lifetable
lx = adjust(s,x,w,radix)$lx
index = which(lx<radix/2)
ifelse(length(index)==0, w, index[1]+x-1)
}
mode <- function(x, w=global$w, radix=global$radix, s=0) {
which.max(adjust(s,x,w,radix)$dx)+x
}
sdtx <- function(x, w=global$w) { # integrate survival function
sqrt(2*integrate(function(t) t*tpx(t,x), 0, w-x)$value -
integrate(tpx, 0, (w-x))$value^2)
}
sdkx <- function(x,w=global$w, radix=global$radix) { # discrete case
lx = createLifeTable(x,w,radix)$lx
k = 2:(w-x)
sqrt(sum(lx[k]*(2*k - 3))/lx[1] -
(sum(lx[k])/lx[1])^2)
}
adjust <- function(s, x=global$x, w=global$w, radix=global$radix) { # reduce lt by s years
if(s==0) return(createLifeTable(x,w,radix))
lt = createLifeTable(x,w,radix)
lt$qx = lt$qx * Rx(x:(w-1), w)^(s+x:(w-1)-x)
lt$px = 1 - lt$qx
lt$dx = radix*c(1, cumprod(lt$px))[1:(w-x)]*lt$qx
lt$lx = radix-c(1,cumsum(lt$dx))[1:(w-x)]
lt
}
lt = createLifeTable(0) # newborns
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
ltx = createLifeTable(global$x) # aged 25
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
ltx10 = adjust(10) # 10-year reduction
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
ltx25 = adjust(25) # 25-year reduction
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
par(mfrow=c(2,2))
plot(tpx,0,(global$w-global$x),ylab="tpx",xlab="T(x)",main="Survival Function")
plot(seq(0:(global$w-global$x-1)), ltx$dx, type='l', main="Curve of Deaths", xlab="T(x)", ylab="Deaths")
output = data.frame(quantities=c("ncurtate(n)", "ncomplete(n)", "median(x)", "mode(x)",
"sdtx(x)", "sdkx(x)", "median(x, s=10)", "mode(x, s=10)", "median(x, s=25)",
"mode(x, s=25)"), values=c(ncurtate(global$n), ncomplete(global$n), median(global$x),
mode(global$x), sdtx(global$x), sdkx(global$x), median(global$x, s=10),
mode(global$x, s=10), median(global$x, s=25), mode(global$x, s=25)))
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
## Warning in tpx(seq(0, w - x - 2), x) * kdeferredqx(1, seq(x, w - 1)):
## longer object length is not a multiple of shorter object length
print(output) # print(lt); print(ltx); print(ltx10); print(ltx25)
## quantities values
## 1 ncurtate(n) 27.41683
## 2 ncomplete(n) 27.71630
## 3 median(x) 58.00000
## 4 mode(x) 26.00000
## 5 sdtx(x) 20.87633
## 6 sdkx(x) 20.87555
## 7 median(x, s=10) 60.00000
## 8 mode(x, s=10) 100.00000
## 9 median(x, s=25) 60.00000
## 10 mode(x, s=25) 100.00000