#__________________________________________________
#test file
require(lifecontingencies)
#__________________________________________________
#vectorize function for testing purposes
pXt <- Vectorize(pxt, "x")
qXt <- Vectorize(qxt, "x")
pxT <- Vectorize(pxt, "t")
qxT <- Vectorize(qxt, "t")
tol <- (.Machine$double.eps)^(1/3)
#see file R/testfunc-demographic.R
getpxqx <- lifecontingencies:::getpxqx
gettpx <- lifecontingencies:::gettpx
gettqx <- lifecontingencies:::gettqx
getkpxqxk <- lifecontingencies:::getkpxqxk
#US life table 1979-1981 p 60 Actuarial mathematics
x <- 0:109
lx <- c(100000,
98000+ c(740, 648, 584, 535, 495, 459, 426, 396, 370, 347, 328, 309, 285, 248, 196, 129, 047),
97000+ c(953, 851, 741, 623, 499, 370, 240, 110),
96000+ c(982, 856, 730, 604, 477, 350, 220, 88),
95000+ c(951, 808, 655, 492, 317, 129),
94000+ c(926, 706, 465, 201),
93000+ c(913, 599, 256),
92000+ c(882, 472, 21), 91526, 90986, 90402, 89771, 89087, 88348, 87551, 86695, 85776, 84789,
83726, 82726, 81348, 80024, 78609, 77107, 75520, 73846, 72082, 70218, 68248, 66165, 63972, 61673,
59279, 56799, 54239, 51599, 48878, 46071, 43180, 40208, 37172, 34095, 31012, 27960, 24961, 22038,
19235, 16598, 14154, 11908, 9863, 8032, 6424, 5043, 3884, 2939, 2185, 1598, 1150, 815, 570, 393,
267, 179, 119, 78, 51, 33)
usLT7981 <- data.frame(x=x, lx=lx)
USLT7981 <- new("lifetable", x=x, lx=lx)
#__________________________________________________
#living prob (one-year)
stopifnot(
sum(abs(
pXt(USLT7981, x=head(x, -1), t=1) - getpxqx(usLT7981$lx)[, "px"])
) < tol
)
#living prob (t-year) at age x
stopifnot(
all(
sapply(c(5, 10, 15, 20, 25, 30, 60, 70), function(k)
sum(abs(pxT(USLT7981, x=k, t=0:(max(x)-k)) - getkpxqxk(usLT7981$lx, x=k)[, "kpx"]))
) < tol)
)
#fractional year living prob
myt <- 1:30/8
stopifnot(
sum(abs(
cbind(
pxT(USLT7981, x=10, t=myt, fractional="linear") - gettpx(usLT7981$lx, x=10, k=myt, fractional="linear"),
pxT(USLT7981, x=10, t=myt, fractional="hyperbolic") - gettpx(usLT7981$lx, x=10, k=myt, fractional="balducci"),
pxT(USLT7981, x=10, t=myt, fractional="constant force") - gettpx(usLT7981$lx, x=10, k=myt, fractional="constant")
)
)) < tol
)
#fractional year death prob
stopifnot(
sum(abs(
cbind(
qxT(USLT7981, x=10, t=myt, fractional="linear") - gettqx(usLT7981$lx, x=10, k=myt, fractional="linear"),
qxT(USLT7981, x=10, t=myt, fractional="hyperbolic") - gettqx(usLT7981$lx, x=10, k=myt, fractional="balducci"),
qxT(USLT7981, x=10, t=myt, fractional="constant force") - gettqx(usLT7981$lx, x=10, k=myt, fractional="constant")
)
)) < tol
)
#__________________________________________________
#death prob (one-year)
stopifnot(
sum(abs(qXt(USLT7981, x=head(x, -1), t=1) - getpxqx(usLT7981$lx)[, "qx"])) < tol
)
#death prob (one-year) at age x+k
stopifnot(
sum(
sapply(c(5, 10, 15, 20, 25, 30), function(k)
sum(abs(qXt(USLT7981, x=k+0:(max(x)-k), t=1) - getkpxqxk(usLT7981$lx, x=k)[, "qxplusk"]))
)) < tol
)
#__________________________________________________
#residual life expectancy
ex <- function(x)
sum(getkpxqxk(usLT7981$lx, x=x)[-1, "kpx"])
stopifnot(
sum(
sapply(0:20, function(x) abs(ex(0) - exn(USLT7981, 0)))
) < tol
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.