\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{STA 2101/442 Assignment Five}}\footnote{Copyright information is at the end of the last page.}
\vspace{1 mm}
\end{center}

\noindent
Please bring your R printouts to the quiz. The non-computer questions are just practice for the quiz, and are not to be handed in, though you may use R as a calculator. Bring a real calculator to the quiz.

\begin{enumerate} 

%%%%%%%%%% Likelihood 1 %%%%%%%%%%%%

\item \label{numle} For each of the following distributions and associated data sets, obtain the maximum likelihood estimate numerically with \texttt{R}. \emph{Bring your printout for each problem to the quiz}.; you may be asked to hand it in. There are links to the data from the course web page in case the ones from this document do not work. 
    \begin{enumerate}

        \item $f(x) = \frac{1}{\pi[1+(x-\theta)^2]}$ for $x$ real, where $-\infty < \theta < \infty$. Data: 
                % 50% mixture of Cauchy(-5) and Cauchy(-5): Two local maxima
                % +4.263357 and -3.719397, global at latter. Signs of data switched from 2011, 
                % which should make it more interesting.

\begin{verbatim}        
-3.77  -3.57  4.10  4.87  -4.18  -4.59  -5.27  -8.33  5.55  -4.35  -0.55  
5.57 -34.78  5.05  2.18  4.12  -3.24  3.78  -3.57  4.86
\end{verbatim}        
        
You might want to read the data from 
\href{http://www.utstat.toronto.edu/~brunner/appliedf13/code_n_data/hw/cauchy.data}
     {\texttt{cauchy.data}}. For this one, try at least two different starting values and \emph{plot the minus log likelihood function!}

%

        \item \label{beta} $f(x) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)} x^{\alpha-1}(1-x)^{\beta-1}$ for $0<x<1$, where $\alpha>0$ and $\beta>0$.  Data: % Beta(10,20) Warn bounds 
        
\begin{verbatim}        
0.45 0.42 0.38 0.26 0.43 0.24 0.32 0.50 0.44 0.29 0.45 0.29 0.29 0.32 0.30
0.32 0.30 0.38 0.43 0.35 0.32 0.33 0.29 0.20 0.46 0.31 0.35 0.27 0.29 0.46
0.43 0.37 0.32 0.28 0.20 0.26 0.39 0.35 0.35 0.24 0.36 0.28 0.32 0.23 0.25
0.43 0.30 0.43 0.33 0.37
\end{verbatim}        
        
You might want to read the data from 
\href{http://www.utstat.toronto.edu/~brunner/appliedf13/code_n_data/hw/beta.data}
     {\texttt{beta.data}}.
If you are getting a lot of warnings, maybe it's because the numerical search is leaving the parameter space. If so, try \texttt{help(nlminb)}.

        \item $f(x) = \frac{\theta e^{\theta(x-\mu)}}{(1+e^{\theta(x-\mu)})^2}$ for $x$ real, where $-\infty < \mu < \infty$ and $\theta > 0$.  Data: 
\begin{verbatim}        
4.82  3.66  4.39  1.66  3.80  4.69  1.73  4.50  9.29  4.05  4.50 -0.64  1.40
4.18  2.70  5.65  5.47  0.55  4.64  1.19  2.28  7.16  4.80  3.19  2.33  2.57
2.31  0.35  2.81  2.35  2.52  3.44  2.71 -1.43  7.61  0.93  2.52  6.86  6.14
4.37  3.79  5.04  4.50  1.92  3.25 -0.06  2.81  3.09  2.95  3.69
\end{verbatim}  
% Logistic mu = 3, theta=1, median 3.22, mle = 
% $estimate
% [1] 3.3392004 0.8760226      
You might want to read the data from 
\href{http://www.utstat.toronto.edu/~brunner/appliedf13/code_n_data/hw/mystery.data}
     {\texttt{mystery.data}}. \vspace{2mm}

\newpage

        \item $f(x) = \frac{1}{m!} e^{-x} x^m$ for $x>0$, where the unknown parameter $m$ is a positive integer.  \emph{This means your estimate will be an integer.} Data: 
% Gamma m=5 so alpha=6,beta=1

\begin{verbatim}        
 8.34  7.65  6.72  3.84  7.12  1.88  5.07  2.69  4.50  5.78  4.88  5.23  6.17
11.76  7.84  5.87  5.23  6.55  8.34  5.35  4.98 13.81  8.62  7.88  6.34  5.16
 6.64  4.35  6.77  5.83  5.85  2.46  8.33  3.74  5.10  3.95  7.84  4.70  6.09
 5.23  1.44  6.11  4.88  7.24  7.89  8.98  1.78  5.46  5.34  4.25
\end{verbatim}        
        
You might want to read the data from 
\href{http://www.utstat.toronto.edu/~brunner/appliedf13/code_n_data/hw/gamma.data}
     {\texttt{gamma.data}}.
        
    \end{enumerate}
    
For each distribution, be able to state (briefly) why differentiating the log likelihood and setting the derivative to zero does not work. For the computer part, bring to the quiz one sheet of printed output for each distribution. The sheets should be separate, because you may hand only one of them in.
Each printed page should show the following, \emph{in this order}. 
\begin{itemize}
    \item Definition of the function that computes the likelihood, or log likelihood, or minus log likelihood or whatever.
    \item How you got the data into R -- probably a \texttt{scan} statement.
    \item Listing of the data for the problem. 
    \item The \texttt{nlm} or \texttt{nlminb} statement and resulting output.
    \item For the Cauchy example, a plot of the minus log likelihood.
\end{itemize}

\item \label{lrtest} For the data of Problem~\ref{beta}, conduct a large-sample likelihood ratio test of $H_0:\alpha=\beta$, using \texttt{R}. Your printout should display the value of $G^2$, the degrees of freedom and the $p$-value. Do you reject $H_0$ at the $0.05$ significance level? If yes, which parameter seems to be higher based on $\widehat{\alpha}$ and $\widehat{\beta}$? \emph{Bring your printout to the quiz}.

\item \label{LRmultinomial} For a random sample from a Multinomial$(1,\boldsymbol{\theta})$ distribution, show that the likelihood ratio test statistic can be written as
\begin{displaymath}
    G^2 = 2n \sum_{j=1}^c \overline{Y}_j\log 
    \left(\frac{\overline{Y}_j}  {\widehat{\theta}_j}\right),
\end{displaymath} 
where $\widehat{\theta}_j$ is the \emph{restricted} maximum likelihood estimate of $\theta_j$. That is, it's the MLE under $H_0$. You may use without proof the fact that the unrestricted MLE is $\overline{Y}_j$.

\newpage

    \item \label{soap} In yet another marketing study, consumers are given samples of three different laundry detergents, which we will call ``A," ``B" and ``C." The labels that the people see are randomly assigned. Suppose $n=150$ consumers participate in the product test. Thirty-eight prefer ``A", 55 prefer ``B,", and 57 prefer ``C." 
    \begin{enumerate}
        \item Test the null hypothesis of equal preference using a large-sample likelihood ratio test. Of course you may use Question~\ref{LRmultinomial}.
        \item Test the same null hypothesis using a Wald test. I used numerical maximum likelihood even though I knew the MLE.
    \end{enumerate}
In both cases, use \texttt{R}. Your printout should display the value of the test statistic, the degrees of freedom and the $p$-value. Do you reject $H_0$ at the $0.05$ significance level? \emph{Bring your printout to the quiz}.

\item \label{mvn} Based on a random sample of size $n$ from a $p$-dimensional multivariate normal distribution, derive a formula for the large-sample likelihood ratio test statistic $G^2$ for the null hypothesis that $\boldsymbol{\Sigma}$ is diagonal (all covariances between variables are zero). You may use material on the slides from the multivariate normal lecture (including a useful version of the likelihood function), which will be provided with the quiz if this question is asked. You may also use without proof the fact that the unrestricted MLE is 
$\widehat{\theta}=(\overline{\mathbf{x}},\boldsymbol{\widehat{\Sigma}})$,
where 
        \begin{displaymath}
            \boldsymbol{\widehat{\Sigma}} = 
            \frac{1}{n}\sum_{i=1}^n (\mathbf{x}_i-\overline{\mathbf{x}}) 
                        (\mathbf{x}_i-\overline{\mathbf{x}})^\prime .
        \end{displaymath}
Hint: Because zero covariance implies independence for the multivariate normal, the joint density is a product of marginals under $H_0$. % From 431s09 Assignment 4, which has a lot more detail.

\item This question uses the SAT data again:  
\href{http://www.utstat.toronto.edu/~brunner/appliedf13/code_n_data/hw/sat.data}
     {\texttt{sat.data}}. There is a link on the course web page in case the one in this document does not work.

Using R, carry out a large-sample likelihood ratio test to determine whether there are any non-zero covariances among the three variables; guided by the usual $\alpha=0.05$ significance level, what do you conclude?  Are the three variables all independent of one another? Answer Yes or No.
\emph{Bring your printout to the quiz}.




\end{enumerate}




% \vspace{25mm}

\noindent
\begin{center}\begin{tabular}{l}
\hspace{6.5in} \\ \hline
\end{tabular}\end{center}
This assignment was prepared by  \href{http://www.utstat.toronto.edu/~brunner}{Jerry Brunner},
Department of Statistics, 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:
\href{http://www.utstat.toronto.edu/~brunner/oldclass/appliedf13} {\texttt{http://www.utstat.toronto.edu/$^\sim$brunner/oldclass/appliedf13}}

\end{document}

# Beta work

rm(list=ls())
x <- scan("http://www.utstat.toronto.edu/~brunner/appliedf13/code_n_data/hw/beta.data") 

bll = function(ab,datta) # - Loglike of beta
    { bll = -sum(dbeta(datta,ab[1],ab[2],log=T)); bll }

# nlm(bll,c(1,1),datta=x) # Works, with warnings.

fit1 <- nlminb(c(1,1),objective=bll,lower=c(0,0),datta=x); fit1

bll0 = function(theta,datta) # - Loglike of beta under H0: alpha=beta
    { bll0 = -sum(dbeta(datta,theta,theta,log=T)); bll0 }

fit0 = nlminb(1,objective=bll0,lower=0,datta=x); fit0

# nlm(bll0,p=1,datta=x) # Works with warnings

# G^2 is twice difference between (minus) log likelihoods)
G2 = fit0$objective-fit1$objective; G2
df=1
pval = 1-pchisq(G2,df); pval






# Old

BLL <- 
     {
     n <- length(datta)
     BLL <- n*lgamma(ab[1]) + n*lgamma(ab[2]) - n*lgamma(sum(ab)) -
     (ab[1]-1)*sum(log(datta)) - (ab[2]-1)*sum(log(1-datta))
     if(ab[1] <= 0) BLL <- Inf ;  if(ab[2] <= 0) BLL <- Inf
     BLL
     }


# Output

> rm(list=ls())
> x <- scan("http://www.utstat.toronto.edu/~brunner/appliedf13/code_n_data/hw/beta.data") 
Read 50 items
> 
> bll = function(ab,datta) # - Loglike of beta
+     { bll = -sum(dbeta(datta,ab[1],ab[2],log=T)); bll }
> 
> # nlm(bll,c(1,1),datta=x) # Works, with warnings.
> 
> fit1 <- nlminb(c(1,1),objective=bll,lower=c(0,0),datta=x); fit1
$par
[1] 13.96757 27.27780

$objective
[1] -60.26451

$convergence
[1] 0

$iterations
[1] 19

$evaluations
function gradient 
      20       46 

$message
[1] "relative convergence (4)"

> 
> bll0 = function(theta,datta) # - Loglike of beta under H0: alpha=beta
+     { bll0 = -sum(dbeta(datta,theta,theta,log=T)); bll0 }
> 
> fit0 = nlminb(1,objective=bll0,lower=0,datta=x); fit0
$par
[1] 3.796628

$objective
[1] -18.13277

$convergence
[1] 0

$iterations
[1] 9

$evaluations
function gradient 
      10       12 

$message
[1] "relative convergence (4)"

> 
> # nlm(bll0,p=1,datta=x) # Works with warnings
> 
> # G^2 is twice difference between (minus) log likelihoods)
> G2 = fit0$objective-fit1$objective; G2
[1] 42.13173
> 
> df=1
> pval = 1-pchisq(G2,df); pval
[1] 8.532719e-11




# Try that multinomial 38 55 57 numerically
mllm = function(theta)
    { n1 = 38; n2 = 55; n3 = 57 
    mllm = - n1 * log(theta[1]) -  n2 * log(theta[2]) - 
             n3 * log(1-theta[1]-theta[2])
    mllm
    }
ybar = c(38/150,55/150); ybar

nlm(mllm,ybar,hessian=T)
# Maybe you get warnings with other starting values, but you still come to y-bar