#
#  本題列出2009年上學期 R 的修勻程式檔案
#
###############################################################################
#
# (1) MWA graduation: 
#     --> 輸入「死亡率」、「修勻係數」(由參數 n 及 z 決定！)
#
mwa=function(data,coef) {
 n=nrow(data)
  m=(length(coef)-1) 
   ux=data[,2]
    vx=rep(0,n) 
     for (i in (m+1):(n-m)) { 
      for (j in (-m:m)) { 
       vx[i]=vx[i]+coef[abs(j)+1]*ux[i+j]
       }
      }
     k=length(vx) 
   vx[1:m]=NA 
  vx[(k-m+1):k]=NA    
 return(vx)
}
# 
#  以下為四種常用的係數，分別為n=2, 4 與 z=0, 3 四種交叉
coef1=c(.4857,.3429,-.0857)
coef2=c(.5594,.2937,-.0734)
coef1=c(.2554,.2338,.1688,.0606,-.0909)
coef2=c(.3311,.2666,.1185,-.0099,-.0407)
#
#  使用方法：資料(data)輸入需為兩行，第二行儲存ux；
#            修勻係數(coef)輸入上述之類的係數，但注意損失觀察值的個數
#
data=matrix(scan("g:\\graduate951(assg1)5.txt"),ncol=7,byrow=T)
age=data[,1]
 ux=data[,3]/data[,2]
  v1=mwa(cbind(age,ux),coef1)
   v2=mwa(cbind(age,ux),coef2)
    v3=mwa(cbind(age,ux),coef3)
     v4=mwa(cbind(age,ux),coef4)
      xx=cbind(age,age,age,age,age)
     yy=cbind(ux,v1,v2,v3,v4)
   matplot(xx,yy,type="l",col=c(4,1,2,1,2),lty=c(3,1,1,4,4),xlab="Age",ylab="qx",main="第一次作業第五題(MWA修勻)")
  legend(65,0.05,c("Raw","MWA(n=2,z=0)","MWA(n=2,z=3)","MWA(n=4,z=0)","MWA(n=4,z=3)"),col=c(4,1,2,1,2),lty=c(3,1,1,4,4))
#
#
###############################################################################
#
# (2) Whittaker graduation:
#     --> 輸入「年齡」「暴露數」「死亡人數」共三行死亡率資料、「加權數」、「參數z」、「參數h」  
#
wh=function(data,W,z,h=0.1)
{
  n=nrow(data)	# the length of the data
  age=data[,1]
  ux=data[,2]
  K=matrix(rep(0,n*(n-z)),n-z,n)
  for (j  in  1: (n-z))  {
    for ( i  in  0:z)  {
     K[j,i+j]=(-1)^(z-i)*gamma(z+1)/(gamma(z-i+1)*gamma(i+1))  
      } 
     }
  vx=solve((W+h*crossprod(K,K)))%*%W%*%ux
  return(vx)
}  
#
#  使用方法：資料(data)輸入需為兩行，第二行儲存ux，第一行為年齡；
#            W 矩陣為加權矩陣，屬於對角矩陣，建議代入各年齡人數(或暴露數)；
#            修勻參數(z, h)的選擇建議為 z=3, h為平均觀察人數(或平均暴露數)
#  註：以下的「data」儲存「年齡」「暴露數」「死亡人數」共三行死亡率資料
# 
age=data[,1]
 ux=data[,3]/data[,2]
  tt=cbind(age,ux)
  W=diag(data[,2])
   v1=wh(tt,W,3,h=10)
    v2=wh(tt,W,3,h=100)
     v3=wh(tt,W,3,h=1000)
      v4=wh(tt,W,3,h=10000)
      xx=cbind(age,age,age,age,age)
     yy=cbind(ux,v1,v2,v3,v4)
   matplot(xx,yy,type="l",col=c(4,1,2,1,2),lty=c(3,1,1,4,4),xlab="Age",ylab="qx",main="第一次作業第五題(Whittaker修勻, z=3)")
  legend(65,0.05,c("Raw","Whittaker, h=10","Whittaker, h=100","Whittaker, h=1000","Whittaker, h=10000"),col=c(4,1,2,1,2),lty=c(3,1,1,4,4))
# --> Choose h=10000
#