% Sample Question document for STA312
\documentclass[12pt]{article} 
%\usepackage{amsbsy} % for \boldsymbol and \pmb 
%\usepackage{graphicx} % To include pdf files!
\usepackage{amsmath}
\usepackage{amsbsy}
\usepackage{amsfonts}
\usepackage[colorlinks=true, pdfstartview=FitV, linkcolor=blue, citecolor=blue, urlcolor=blue]{hyperref} % For links
\usepackage{fullpage}
%\pagestyle{empty} % No page numbers


\begin{document}
%\enlargethispage*{1000 pt} 

\begin{center}   
{\Large \textbf{Sample Questions: Survival and Hazard Functions}}%\footnote{}
\vspace{1 mm}

STA312 Fall 2023. Copyright information is at the end of the last page.
\end{center}

\vspace{5mm}

\noindent
For all these questions, $T$ is a continuous random variable with $P(T>0)=1$, density $f(t)$ and cumulatve distribution function $F(t) = P(T \leq t)$. 

\begin{enumerate} 

\item The survival function is $S(t) = P(T>t)$. Prove 
{ \Large $\displaystyle E(T) = \int_0^\infty S(t) \, dt$. }

\pagebreak %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\item The hazard function is defined by 
{ \Large $\displaystyle  h(t) =  \lim_{\Delta \rightarrow 0} \frac{P(t < T < t + \Delta|T>t)}{\Delta}$, }
where $\Delta>0$. Prove
{ \Large $ h(t) = \frac{f(t)}{S(t)}$. }

\pagebreak %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\item Prove {\Large  $\displaystyle S(t) = e^{-\int_0^t h(x) \, dx}$. }

\pagebreak %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\item  {\large Let $T \sim \exp(\lambda)$. Find the hazard function $h(t)$ for $t>0$. }

\pagebreak %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\item Let $T$ have the Pareto density 
{ \large
$f(t|\theta) = \left\{ \begin{array}{ll} % ll means left left
   \frac{\theta}{t^{\theta+1}} & \mbox{for $ t \geq 1$}  \\
   0     & \mbox{for } t<1 
   \end{array}  \right. $ % 
} % End size
    \begin{enumerate}
        \item Find the hazard function $h(t)$ for $t>1$. % theta/t. 

\pagebreak %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        \item Earlier, we found the MLE { \Large $\widehat{\theta}_n = \frac{n}{\sum_{i=1}^n \log t_i}$,} and {\Large $\widehat{v}_n = \frac{\widehat{\theta}_n^2}{n}$. } \vspace{1mm}

        \begin{enumerate}
            \item Give $\widehat{h(t)}$, the maximum likelihood estimate of the hazard function evaluated at a particular time $t>1$. Your answer is a formula involving $t$ and $\widehat{\theta}_n$. \vspace{40mm}
            \item We want a confidence interval for $h(t)$, the hazard function evaluated at a particular time $t>1$. Give formulas for the lower and upper 95\% confidence limits. Show your work. 
        \end{enumerate}



    \end{enumerate}

\pagebreak %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\item Let $T$ have a gamma distribution with parameters $\alpha>0$ and $\lambda>0$. 
        \begin{enumerate}
            \item What is the hazard function? \vspace{50mm}
            \item Using R, plot the hazard function for several values of $\alpha$ and $\lambda$. How do the parameter values influence the shape of the hazard function? 
        \end{enumerate}

\end{enumerate} % End of all the questions

\vspace{100mm}

\noindent
\begin{center}\begin{tabular}{l}
\hspace{6in} \\ \hline
\end{tabular}\end{center}
This assignment was prepared by  \href{http://www.utstat.toronto.edu/~brunner}{Jerry Brunner},
Department of Mathematical and Computational Sciences, University of Toronto. It is licensed under a 
\href{http://creativecommons.org/licenses/by-sa/3.0/deed.en_US}
     {Creative Commons Attribution - ShareAlike 3.0 Unported License}. Use any part of it as you like and share the result freely. The \LaTeX~source code is available from the course website: 
     
\begin{center}
\href{http://www.utstat.toronto.edu/brunner/oldclass/312f23} {\small\texttt{http://www.utstat.toronto.edu/brunner/oldclass/312f23}}
\end{center}

\end{document}

##########################################################################
# Hazard function of Gamma(alpha,lambda) distribution: : h(t) = f(t)/S(t)

# You might think the following would work, but watch.
# Remember that a gamma distribution with alpha=1 is exponential.
# The hazard function of an exponential distribution is a constant lambda.

Time = 1:10

alpha=1; lambda=1
dgamma(Time,shape=alpha,rate=lambda) / (1-pgamma(Time,shape=alpha,rate=lambda))

alpha=1; lambda=2
dgamma(Time,shape=alpha,rate=lambda) / (1-pgamma(Time,shape=alpha,rate=lambda))

alpha=1; lambda=3
dgamma(Time,shape=alpha,rate=lambda) / (1-pgamma(Time,shape=alpha,rate=lambda))

alpha=1; lambda=4
dgamma(Time,shape=alpha,rate=lambda) / (1-pgamma(Time,shape=alpha,rate=lambda))

alpha=1; lambda=10
dgamma(Time,shape=alpha,rate=lambda) / (1-pgamma(Time,shape=alpha,rate=lambda))

# Use logs to get around division by a number close to zero
alpha=1; lambda=10
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
exp(logh)


# Now plot
alpha=1/2;lambda=1
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell
# Notice explicit limits on y to keep all plots on the same scale.

alpha=1/2;lambda=2
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

alpha=1/2;lambda=3
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

# Try lambda small
alpha=1/2;lambda=1/2
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

alpha=1/2;lambda=1/4
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

# Lambda appears to control the overall level of the function, with larger numbers higher.
# Investigate alpha, initially holding lambda fixed at 2

alpha=2;lambda=2
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

alpha=3;lambda=2
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

alpha=4;lambda=2
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

# Increase lambda to 5

alpha=4;lambda=5
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

alpha=3;lambda=5
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

alpha=2;lambda=5
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

# So far, we know that 
    # If alpha=1, h(t) is constant at lambda
    # If alpha > 1, h(t) is increasing and bounded above by lambda
    # If alpha < 1, h(t) is decreasing and bounded below by lambda
    # alpha values close to one may make the big change occur closer to zero.

# Now alpha < 1 again
alpha=1/2;lambda=2
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

alpha=1/10;lambda=2
Time = seq(from=0,to=5,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell

# Is the hazard function asymptotic to lambda?

alpha=1/10;lambda=2
Time = seq(from=0,to=50,length=101)
logh = dgamma(Time,shape=alpha,rate=lambda,log=TRUE) - 
    pgamma(Time,shape=alpha,rate=lambda, lower.tail=FALSE,log.p=TRUE)
Hazard = exp(logh)
tstring = paste('alpha =',alpha,'and lambda =',lambda)
plot(Time,Hazard,type='l',main=tstring,ylim=c(0,5)) # That's a lower case ell
# Add the line y = lambda = 2
xx = c(0,50); yy =c(2,2); lines(xx,yy,lty=2)

# This is no proof, but I am convinced that the hazard function is asymptotic to lambda.