# ("current" or "period") lifetable calculations

ds=read.table("dths_pop_can2001_for_R.txt",header=T) ; ds
attach(ds)

age=1:101

rate.male=d.male/(p.male*1) ; rate.male ; plot(age,rate.male)
 # note that each rate is an incidence density, since
 # the denominator = (mid-year pop) x 1-year  
 
p.male = exp(- rate.male * 1) ; p.male ; plot(age,p.male)
 # rate.male * 1 is the integral. 
 # p is the surviving fraction, of those alive at start of year
 # this is not the classical way to calculate p
 
q.male = 1 - p.male ; q.male ; plot(age,q.male) 

l.male = round(100000*cumprod(p.male)) ; 
  l.male ; plot(age,l.male) ;
 # the number living at each birthday
 # remember one can only die once, but one is counted as a
 # a survivor at each birthday.. thus 'fraction of a fraction ..'
 
L.male = round( ( c(100000,l.male[1:100]) + l.male) / 2 ) ; L
  sum(L) / 100000
 # approx no. of py lived in this year
 
d.male = c(100000,l.male[1:100]) - l; 
  d.male; sum(d.male) ; plot(age,d.male)
  sum((age-0.5)*d.male)/sum(d.male)
 # no. of deaths in this year
 
T.male = (cumsum(L[101:1]))[101:1] ; T.male
 #  remaining person years - cumulates from 101 to 1;
 #  then reverses the order to go forward again
 
e.male = T.male / c(100000,l.male[1:100]) ; 
   e.male; plot(age,e.male) ;
 # average (expectation) of remaining person years