actuarial <- function(time,status,strata=NULL,weight=NULL,intervals){
   if(is.null(weight)) weight<-rep(1,length(time))
   if(is.null(strata)) strata<-rep(1,length(time))
   newstrat<-match(strata,unique(strata))
   nstrat<-length(unique(strata))
   nint<-length(intervals)-1
   n.event<-n.cens<-n.risk<-array(NA,c(nint+1,nstrat))
   for(j in seq(length(unique(strata)))){
      for(i in seq(nint+1)){
         f<-if(i<=nint) (time<=intervals[i+1]) else rep(1,length(time))
         n.risk[i,j]<-sum(weight*(time>intervals[i])*(newstrat==j))
         n.event[i,j]<-sum(weight*(time>intervals[i])*(newstrat==j)*status*f)
         n.cens[i,j]<-sum(weight*(time>intervals[i])*(newstrat==j)*(1-status)*f)
      }
   }
   n.risk.adj<-n.risk-n.cens/2
   p<-1-n.event/n.risk.adj
   p[n.risk.adj==0]<-1
   q<-1-p
   survival<-exp(apply(log(rbind(1,p[-nint-1,,drop=F])),2,cumsum))
   dens<--apply(rbind(survival,NA),2,diff)/outer(c(diff(intervals),NA),
      rep(1,nstrat))
   haz<-2*dens/(rbind(survival[-1,,drop=F],NA)+survival)
   vc<-rbind(0,(q/(n.risk.adj*p))[-nint-1,,drop=F])
   vc[n.risk.adj==0]<-0
   lsds<-sqrt(apply(vc,2,cumsum))
   sds<-survival*lsds
   return(as.data.frame(list(
      stratum=unique(strata)[rep(seq(nstrat),rep(nint+1,nstrat))],
      time=as.vector(outer(c((intervals[-1]+intervals[-(nint+1)])/2,intervals[nint+1]+1),
         rep(1,nstrat))),
      n.risk=as.vector(n.risk),
      n.risk.adj=as.vector(n.risk.adj),n.event=as.vector(n.event),
      survival=as.vector(survival),dens=as.vector(dens),haz=as.vector(haz),
      se.survival=as.vector(sds)
   )))
}