# R/DiscreteHazRateEst.R In NPHazardRate: Nonparametric Hazard Rate Estimation

#### Documented in baseCparamCalculationlambdahatnsfpower.matrixSemiparamEstSmoothedEstimateTm

lambdahat<-function(xin, cens, xout)
#non smooth estimate of the discrete hazard function
{
n<-length(xout)
mat<-sapply(1:n, function(i, xout, xin, cens)
{
xoutuse<-xout[i]
whichtk<-which(xoutuse == xin) ## index of those T_j = t_k
xinuse<-xin[whichtk]
deltause<-cens[whichtk]
num<- length(which(deltause !=	0))	#equivalent to \sum_{i=1}^n I(T_i = t_k)\Delta_i
den<-length(xinuse)
c(num, den)
}, xout, xin, cens)
mat<-matrix(mat, ncol=2, byrow=TRUE)

numerator<-mat[,1]
denuse<-mat[,2]
m <- length(denuse)
denominator<-sapply(1:m,  function(j, denuse, m) sum(denuse[j:m]), denuse=denuse, m=m) #denominator: \sum_{j \geq k} \sum_{k=1}^n I(T_j=t_k)
numerator/(denominator+1)
}

nsf<-function(xin, cens, xout)
{
n<-length(xout)
mat<-sapply(1:n, function(i, xout, xin, cens)
{
xoutuse<-xout[i]
whichtk<-which(xoutuse == xin) ## index of those T_j = t_k
length(xin[whichtk]	)
}, xout, xin, cens)
#mat<-matrix(mat, ncol=2, byrow=T)
m <- length(mat)
denominator<-sapply(1:n,  function(j, mat, n) sum(mat[j:n]), mat=mat, n=n)
denominator
}

"Habbema"<-function(xin, x)
{
n<-length(x)
nsuse<-nsf(xin, xin, x)
#cat(nsuse)
huse<- .8 * 0.08 * (nsuse/length(xin))^1.5
#cat(huse)
mat<-sapply(1:n, function(i, x, xin, huse)
{
xoutuse<-x[i]
hh<-huse[i]
(1-hh)^((xoutuse-xin)^2)
}, x, xin, huse)
mat
}

"TutzPritscher"<-function(xin, cens, xout)
{
mle<-lambdahat(xin, cens, xout)
nsuse<-nsf(xin, cens, xin)
#nsuse<-matrix(nsuse, nrow=length(xin), ncol=length(xout))
huse<-bw.nrd(xin)
kernel<-Habbema(xin, xout)
weightsuse<- nsuse * kernel
weights<-weightsuse/colSums(weightsuse)
colSums(mle * weights, na.rm=TRUE)
}

Tm<-function(tk, xout, distribution, par1, par2)
# This is used to calculate max(t_k; 1-\sum_{l=0}^k \eta_l > \epsilon), \epseilon>0
{
SumDistrTk<-PdfSwitch(xout, distribution, par1, par2)
n<-length(SumDistrTk)
SumTK<-sapply(1:n,  function(i, SumDistrTk, n) sum(SumDistrTk[i:n]), SumDistrTk=SumDistrTk, n=n) #\sum_{l=0}^n \eta_l
Tmh<- 1-SumTK  #1-\sum_{l=0}^n \eta_l

maxTmh<-max(which(Tmh == max(Tmh)) )#Tmh[Tmh>= 0 & Tmh<= min(Tmh[Tmh>=0])] # max{1-\sum_{l=0}^n \eta_l >0}
tk[maxTmh] # find the gridpoint tk which corresponds to Tmh, return the result.
}

CparamCalculation<-function(gamparam, VehHazard)
# returns the C smoothing parameter calculated as
# C= gamma/max_{k \geq 0} ( \lambda(t_{k-1}) + \lambda(t_k) + \lambda(t_{k+1}) )
{
VehHazard<-c(0, VehHazard) # \lambda(t_{-1}) is set to 0 here
n<-length(VehHazard)-2
Lsums<-sapply(1:n, function(i, VehHazard) sum(VehHazard[seq(i,i+2, by=1)]), VehHazard) #calculate \lambda(t_{k-1}) + \lambda(t_k) + \lambda(t_{k+1})

gamparam/max(Lsums) #return C
}

##################################################################################################################
##############alternative DetermineSCprod<-function(NonParEst, VehHazard, Hvector, n, gammapar) function #########

# DetermineSCprod<-function(NonParEst, VehHazard, Hvector, n, gammapar)
#   #this finds SC = \gamma((n+1) \hat B_1)^{-1} \hat V_1
#   # n = number of obs, gammapar = sum of vehicle haz at xout (computed elsewhere)
# {
#   VehHazard<-c(0, VehHazard)
#   NonParEst<-c(0, NonParEst)
#   n<-length(VehHazard)-2
#   B1hatElements<-sapply(1:n, function(i, VehHazard, NonParEst)
#     ((NonParEst[i] + NonParEst[i+2])* VehHazard[i+1] - (VehHazard[i] + VehHazard[i+2])* NonParEst[i+1] )^2
#     , VehHazard, NonParEst) #calculate (\hat \lambda(t_{k-1}) + \hat \lambda(t_{k+1})) \lambda(t_k) - ...)^2
#
#   B1hat<-sum(B1hatElements)
#
#   V1hatElements<-sapply(1:n, function(i, VehHazard, NonParEst)
#     NonParEst[i+1]*(1-NonParEst[i+1]) *(VehHazard[i] + VehHazard[i+2])/Hvector[i]
#     , VehHazard, NonParEst)
#
#   V1hat<- sum(V1hatElements	)
#   gammapar *  V1hat / ((n+1)*B1hat)	#estimated optimal product for S and C
#
# }

#######################################################################################################################

# DetermineSCprod<-function(NonParEst, VehHazard, Hvector, nn, gammapar)
#   #this finds SC = \gamma((n+1) \hat B_1)^{-1} \hat V_1
#   # n = number of obs, gammapar = sum of vehicle haz at xout (computed elsewhere)
# {
#   VehHazard<-c(0, VehHazard)
#   NonParEst<-c(0, NonParEst)
#   n<-length(VehHazard)-2
#   B1hatElements<-sapply(1:n, function(i, VehHazard, NonParEst)
#     ( sum(NonParEst[seq(i,i+2, by=2)])* VehHazard[(i+i+2)/2] - sum(VehHazard[seq(i,i+2, by=2)])* NonParEst[(i+i+2)/2])^2
#     , VehHazard, NonParEst) #calculate (\hat \lambda(t_{k-1}) + \hat \lambda(t_{k+1})) \lambda(t_k) - ...)^2
#
#   B1hat<-sum(B1hatElements)
#   V1hatElements<-sapply(1:n, function(i, VehHazard, NonParEst, Hvector){
#     NonParEst[(i+i+2)/2]*(1-NonParEst[(i+i+2)/2]) *sum(VehHazard[seq(i,i+2, by=2)])/Hvector[i]}
#     , VehHazard, NonParEst, Hvector)
#   V1hat<- sum(V1hatElements)
#   (gammapar *  V1hat) / ((nn+1)*B1hat)	#estimated optimal product for S and C
# }

power.matrix <- function(M, n){
#Raise matrix M to the n'th power
powers <- base(n, 2)
result <- diag( nrow(M) )
for(i in 1:length(powers)) {
if(powers[i])
result <- result %*% M
M <- M %*% M
}
result
}

base <- function(m, b){
#Express m as powers of b; m may be vector
maxi <- floor(log(max(m))/log(b))
result <- matrix( 0, length(m), maxi+1 )
for(i in 1:(maxi+1)) {
result[,i] <- m %% b
m <- m %/% b
}
result
}

SmoothedEstimate<-function(NonParEst, VehHazard, gammapar, SCproduct, Cpar)
{
n<-length(VehHazard)-2 # maybe length(VehHazard)-2
Sparam<- floor(SCproduct/Cpar) # S = S*C/C - this may need to change
b0<- gammapar - Cpar*VehHazard[2] # b0 = \gamma - C \lambda_1 (counting starts at zero in the manuscript)
bkplus1<- gammapar - Cpar * head(tail(VehHazard,2),1) # thisis the last value of the Bk vector (not specified in manuscript)
ak <- Cpar * VehHazard
bk <- sapply(1:n, function(i, VehHazard, Cpar, gammapar)
{ gammapar - Cpar * (VehHazard[i] + VehHazard[i+2])}, VehHazard, Cpar, gammapar)
bk<- c(b0,bk, bkplus1)
Wmatrix<-diag(bk)

Wmatrix[ row(Wmatrix) - col(Wmatrix) ==1]<- tail(ak, length(ak)-1)
LambdaMat<- (1/gammapar) * ( diag(1, nrow = length(VehHazard), ncol = length(VehHazard)) %*% Wmatrix)

#################### test #################################
#cat("\n Vehicle Hazard = ", VehHazard, "\n LAMBDA MAT = ", NonParEst %*% power.matrix(LambdaMat,200), "\ndiag= ", diag(power.matrix(LambdaMat,200)),  "\n")
#aaaa<- NonParEst %*% power.matrix(LambdaMat,5)
#plot(1:length(aaaa), VehHazard, type="l")
#lines(1:length(aaaa), aaaa, lty=4)
###########################################################

if(Sparam == 0) LambdaMatrixSpower <- diag(1, nrow=length(VehHazard), ncol= length(VehHazard))
else LambdaMatrixSpower <- power.matrix(LambdaMat, Sparam)

SmoothEst<- matrix(NonParEst, nrow=1, ncol=length(VehHazard)) %*% LambdaMatrixSpower
SmoothEst

}

# VehicleHazardRate<-function(xout, vehicledistr, param1, param2)
# {
#   switch(vehicledistr, binomial=dbinom(xout, param1, param2)/(1-pbinom(xout-1, param1, param2))
#          , geometric = dgeom(xout, param1)/(1-pgeom(xout-1, param1))
#          , poisson = dpois(xout, param1)/(1-ppois(xout-1, param1))
#          , negativebinomial = dnbinom(xout, param1, param2)/(1-pnbinom(xout-1, param1, param2))
#
#   )
# }

SemiparamEst<-function(xin, cens, xout, Xdistr, Udistr, vehicledistr, Xpar1=1, Xpar2=0.5, Upar1=1, Upar2=0.5, vdparam1=1, vdparam2=0.5)
{

TkCutOff<-min(Tm(xout, xout, Xdistr, Xpar1, Xpar2) , Tm(xout, xout, Udistr, Upar1, Upar2), Tm(xout, xout, vehicledistr, vdparam1, vdparam2))
# cat(xout, " ", TkCutOff)

result<-matrix(nrow=length(xout), ncol=2)

Qvector <- 1-cumsum(PdfSwitch(xout-1, Xdistr, Xpar1, Xpar2))
Pvector <- 1-cumsum(PdfSwitch(xout-1, Udistr, Upar1, Upar2))

Hvector <-  Pvector * Qvector
NonparamEst<-lambdahat(xin, cens, xout)

VehicleHazard<-HazardRate(xout, vehicledistr, vdparam1, vdparam2)
GammaParam <- 2 #sum(VehicleHazard)
Cpar <-  0.5# CparamCalculation(GammaParam , VehicleHazard) #upper bound for C 0.5#
SCproduct<- 32* Cpar # DetermineSCprod(NonparamEst, VehicleHazard, Hvector , length(xin), GammaParam)

SmoothEst<- SmoothedEstimate(NonparamEst, VehicleHazard, GammaParam , SCproduct, Cpar )
result[,1]<-xout
result[,2]<-SmoothEst
result
}


## Try the NPHazardRate package in your browser

Any scripts or data that you put into this service are public.

NPHazardRate documentation built on May 2, 2019, 10:24 a.m.