\documentclass{article}\usepackage[]{graphicx}\usepackage[]{xcolor}
% maxwidth is the original width if it is less than linewidth
% otherwise use linewidth (to make sure the graphics do not exceed the margin)
\makeatletter
\def\maxwidth{ %
  \ifdim\Gin@nat@width>\linewidth
    \linewidth
  \else
    \Gin@nat@width
  \fi
}
\makeatother

\definecolor{fgcolor}{rgb}{0.345, 0.345, 0.345}
\newcommand{\hlnum}[1]{\textcolor[rgb]{0.686,0.059,0.569}{#1}}%
\newcommand{\hlstr}[1]{\textcolor[rgb]{0.192,0.494,0.8}{#1}}%
\newcommand{\hlcom}[1]{\textcolor[rgb]{0.678,0.584,0.686}{\textit{#1}}}%
\newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}%
\newcommand{\hlstd}[1]{\textcolor[rgb]{0.345,0.345,0.345}{#1}}%
\newcommand{\hlkwa}[1]{\textcolor[rgb]{0.161,0.373,0.58}{\textbf{#1}}}%
\newcommand{\hlkwb}[1]{\textcolor[rgb]{0.69,0.353,0.396}{#1}}%
\newcommand{\hlkwc}[1]{\textcolor[rgb]{0.333,0.667,0.333}{#1}}%
\newcommand{\hlkwd}[1]{\textcolor[rgb]{0.737,0.353,0.396}{\textbf{#1}}}%
\let\hlipl\hlkwb

\usepackage{framed}
\makeatletter
\newenvironment{kframe}{%
 \def\at@end@of@kframe{}%
 \ifinner\ifhmode%
  \def\at@end@of@kframe{\end{minipage}}%
  \begin{minipage}{\columnwidth}%
 \fi\fi%
 \def\FrameCommand##1{\hskip\@totalleftmargin \hskip-\fboxsep
 \colorbox{shadecolor}{##1}\hskip-\fboxsep
     % There is no \\@totalrightmargin, so:
     \hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}%
 \MakeFramed {\advance\hsize-\width
   \@totalleftmargin\z@ \linewidth\hsize
   \@setminipage}}%
 {\par\unskip\endMakeFramed%
 \at@end@of@kframe}
\makeatother

\definecolor{shadecolor}{rgb}{.97, .97, .97}
\definecolor{messagecolor}{rgb}{0, 0, 0}
\definecolor{warningcolor}{rgb}{1, 0, 1}
\definecolor{errorcolor}{rgb}{1, 0, 0}
\newenvironment{knitrout}{}{} % an empty environment to be redefined in TeX

\usepackage{alltt}
\usepackage[utf8]{inputenc}
\usepackage[spanish]{babel}
\usepackage{geometry}
\usepackage{colonequals}
\usepackage{amsmath}
\usepackage{dsfont}
\usepackage{url}
\newcommand{\apuno}{Halla un estadístico suficiente
  minimal para $\theta=(p,a)$.}
\newcommand{\apdos}{Mediante el método de los
  momentos, halla estimadores para $p$ y $a$ y aplícalos a $\vec x$.}
\newcommand{\aptre}{Halla la estimación
  máximo-verosímil de $p$ y de $a$.  Compara el resultado con el del método de los momentos.}
\newcommand{\apcua}{Calcula mediante \emph{bootstrap}
  un intervalo de confianza para $\mu=p/a$ a nivel $1-\alpha=0.95$.}
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\begin{document}
\raggedright
\begin{enumerate}

\item (4 puntos) Sea $X$ una población con distribución $\gamma(p,a)$
  de la que se obtiene la muestra aleatoria simple $\vec x$ de tamaño
  $20$ siguiente:
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{verbatim}
## c(1.3, 1.2, 3.5, 4, 2.8, 5.1, 1.5, 2.3, 3.8, 4.3, 3, 2, 0.9, 
## 2.5, 1.7, 1.6, 1.1, 1.8, 1.4, 1)
\end{verbatim}
\end{kframe}
\end{knitrout}
Teniendo en cuenta que
$$
f_X(x)=\frac{a^p}{\Gamma(p)}x^{p-1}e^{-ax}I(x>0)\qquad
E(X)=\frac pa\qquad\mbox{Var}(X)=\frac p{a^2}
$$

\begin{enumerate}
\item (1 punto) \apuno
\item (1 punto) \apdos
\item (1 punto) \aptre
\item (1 punto) \apcua
\end{enumerate}
\end{enumerate}
\section*{Resolución}
\begin{enumerate}
\item
  \begin{enumerate}
  \item \apuno

    Para $x>0$, $f(x,p,a)=\frac{a^p}{\Gamma(p)}x^{p-1}e^{-ax}$ $\Longrightarrow$
    $\mathcal L(\vec x,p,a)=\frac{a^{np}}{\Gamma(p)^n}\left(\prod x_i\right)^{p-1}e^{-a\sum x_i}$ $\Longrightarrow$
    $T\colonequals\left(\prod X_i,\sum X_i\right)$ es suficiente por el teorema de factorización, tomando
    $h(\vec x)=1$ y $g\Bigl(\vec t,\vec\theta\Bigr)$ $=$
    $g(t_1,t_2,p,a)$ $=$ $\frac{a^{np}}{\Gamma(p)^n}\left(t_1\right)^{p-1}e^{-a t_2}$, $t_1(\vec x)=\prod x_i$, $t_2(\vec x)=\sum x_i$.
    [Con más rigor, se podría considerar $x$ arbitrario y tomar $h(\vec x) = I(x_{(1)})$, etcétera.]

    Es minimal porque $\frac{\mathcal L(\vec x,p,a)}{\mathcal L(\vec y,p,a)}$ $=$
    $\left(\frac{\prod x_i}{\prod y_i}\right)^{p-1}e^{-a\left(\sum x_i-\sum y_i\right)}$
    es independiente de $(p,a)$ $\iff$ $\frac{\prod x_i}{\prod y_i}=1$ $\land$ $\sum x_i-\sum y_i=0$ $\iff$
    $\left(\prod x_i,\sum x_i\right)=\left(\prod y_i,\sum y_i\right)$ $\iff$ $T(\vec x)=T(\vec y)$.

    Otra forma de verlo sería considerando que la distribución $\gamma$ pertenece a la familia exponencial $2$-paramétrica,
    pues su soporte $(0,\infty)$ no depende de $\theta$ y
    $$f(x,p,a) = \frac{a^p}{\Gamma(p)} \cdot e^{(p-1)\ln x-ax}=c(\theta)h(x)e^{Q_1(\theta)T_1(x)+Q_2(\theta)T_2(x)}$$
    con $c(p,a)=\frac{a^p}{\Gamma(p)}$, $h(x)=1$, $Q_1(p,a)=p-1$, $T_1(x)=\ln x$,  $Q_2(p,a)=-a$ y $T_2(x)=x$, o bien
    $$f(x,\eta_1,\eta_2) = \frac{(-\eta_2)^p}{\Gamma(\eta_1+1)} \cdot e^{\eta_1\ln x+\eta_2x}=c(\eta)h(x)e^{\eta_1T_1(x)+\eta_2T_2(x)}$$
    con la
    parametrización natural $\eta_1=p-1$, $\eta_2=-a$, $\eta=(\eta_1,\eta_2)$, cuyo espacio paramétrico ${\rm H}$ $=$
    $\{(\eta_1,\eta_2)\mid\eta_1>-1,\eta_2<0\}$ contiene un abierto de $\mathds R^2$ (por ejemplo, el abierto $(0,1)\times(-1,1)$ está
    contenido en $\rm H$),
    luego
    $\bigl(\sum T_1(X_i),\sum T_2(X_i)\bigr)$ $=$ $\left(\sum\ln X_i=\ln\prod X_i,\sum X_i\right)$ es suficiente minimal, luego
    también $\left(\prod X_i,\sum X_i\right)$ es suficiente minimal, pues inducen la misma partición ($\ln$ es inyectiva).
  \item \apdos

    $\mu=p/a$, $\sigma^2=p/a^2$ $\Longrightarrow$ $p=\mu^2/\sigma^2$, $a=\mu/\sigma^2$
    $\Longrightarrow$ $\hat p_{\rm MM}=\bar X^2/\hat S^2$, $\hat a_{\rm MM}=\bar X/\hat S^2$
    o, con el estimador sesgado varianza muestral, $\hat p_{\rm MM}=\bar X^2/S^2$, $\hat a_{\rm MM}=\bar X/S^2$
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{mu} \hlkwb{<-} \hlkwd{mean}\hlstd{(x); s2} \hlkwb{<-} \hlkwd{var}\hlstd{(x)}
\hlstd{(pMM} \hlkwb{<-} \hlstd{mu}\hlopt{^}\hlnum{2}\hlopt{/}\hlstd{s2)}
\end{alltt}
\begin{verbatim}
## [1] 3.574151
\end{verbatim}
\begin{alltt}
\hlstd{(aMM} \hlkwb{<-} \hlstd{mu}\hlopt{/}\hlstd{s2)}
\end{alltt}
\begin{verbatim}
## [1] 1.527415
\end{verbatim}
\begin{alltt}
\hlcom{## con el estimador sesgado varianza muestral:}
\hlstd{s2} \hlkwb{<-} \hlkwd{mean}\hlstd{(x}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{-} \hlstd{mu}\hlopt{^}\hlnum{2}
\hlstd{(pMM} \hlkwb{<-} \hlstd{mu}\hlopt{^}\hlnum{2}\hlopt{/}\hlstd{s2)}
\end{alltt}
\begin{verbatim}
## [1] 3.762265
\end{verbatim}
\begin{alltt}
\hlstd{(aMM} \hlkwb{<-} \hlstd{mu}\hlopt{/}\hlstd{s2)}
\end{alltt}
\begin{verbatim}
## [1] 1.607805
\end{verbatim}
\end{kframe}
\end{knitrout}
Con más rigor, se puede partir de los momentos absolutos:
\begin{eqnarray*}
E(X)&=&p/a\implies p=a\cdot E(X)\\
  E(X^2)&=&{\rm Var}(X)+E^2(X)=\frac p{a^2}+\frac{p^2}{a^2}=\frac{p+p^2}{a^2}\\
  &\implies&E(X^2)=\frac{a\cdot E(X)+a^2\cdot E^2(X)}{a^2}\\
  &\implies&a^2\cdot E(X^2)=a\cdot E(X)+a^2\cdot E^2(X)\\
  &\implies&a\cdot E(X^2)=E(X)+a\cdot E^2(X)\\
  &\implies&a=\frac{E(X)}{E(X^2)-E^2(X)}\implies
             \hat a_{\rm MM}=\frac{\bar X}{\overline{X^2}-\bar X^2}=\frac{\bar X}{S^2}\\
    &\implies&p=\frac{E(X)}{E(X^2)-E^2(X)}E(X)\implies
  \hat p_{\rm MM}=\frac{\bar X}{\overline{X^2}-\bar X^2}\bar X=\frac{\bar X^2}{S^2}
\end{eqnarray*}
\item \aptre
  
  $\mathcal L(\vec x,p,a)=\frac{a^{np}}{\Gamma(p)^n}\prod x_i^{p-1}e^{-a\sum x_i}$ $\Longrightarrow$
  $l(a)\colonequals\ln\mathcal L(\vec x,p,a)=np\ln a-n\ln\Gamma(p)+(p-1)\sum\ln x_i-a\sum x_i$ $\Longrightarrow$
  $0=l'(a)=\frac{np}a-\sum x_i$ $\Longrightarrow$
  $ \begin{cases} l''(a)=-\frac{np}{a^2}<0 \text{ (máximo)} \\ \hat a_{\rm MV}=\hat p_{\rm MV}/{\bar X} \end{cases} $

  $l(p)\colonequals\ln\mathcal L(\vec x,p,a)=np\ln a-n\ln\Gamma(p)+(p-1)\sum\ln x_i-a\sum x_i$ $\Longrightarrow$
  $0=l'(p)=n\ln a-n\frac{\Gamma'(p)}{\Gamma(p)}+\sum\ln x_i$
  $=$ $n\ln\frac p{\bar X}-n\psi(p)+\sum\ln x_i$
  $=$ $n\ln p - \ln{\bar X}-n\psi(p)+\sum\ln x_i$
  $=$ $n\ln p -n\psi(p)+\text{constante}$ $\Longrightarrow$
  $\ln p -\psi(p)=\text{constante}$ $\Longrightarrow$ existe a lo sumo una única $p$ que
  lo satisfaga, ya que $\ln$ y $\psi$ son estrictamente
  crecientes\footnote{\url{https://en.wikipedia.org/wiki/Digamma_function}} en
  $(0,\infty)$
  
  Hay que hallar $p$ mediante métodos numéricos, bien como raíz de la ecuación anterior, donde $\psi$=\texttt{digamma},
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{g} \hlkwb{<-} \hlkwa{function} \hlstd{(}\hlkwc{p}\hlstd{) n}\hlopt{*}\hlkwd{log}\hlstd{(p}\hlopt{/}\hlkwd{mean}\hlstd{(x))} \hlopt{-} \hlstd{n}\hlopt{*}\hlkwd{digamma}\hlstd{(p)} \hlopt{+} \hlkwd{sum}\hlstd{(}\hlkwd{log}\hlstd{(x))}
\hlkwd{uniroot} \hlstd{(g,} \hlkwd{c} \hlstd{(}\hlnum{1e-100}\hlstd{,} \hlnum{100}\hlstd{))}
\end{alltt}
\begin{verbatim}
## $root
## [1] 3.997045
## 
## $f.root
## [1] -6.296073e-08
## 
## $iter
## [1] 12
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
\end{verbatim}
\end{kframe}
\end{knitrout}
bien mediante la maximización de la verosimilitud perfilada en $p$:
  $\ln\mathcal L(x,p)=np\ln (p/\bar x)-n\ln\Gamma(p)+(p-1)\sum\ln x_i-np$
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{vp} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{p}\hlstd{)}
    \hlstd{n}\hlopt{*}\hlstd{p}\hlopt{*}\hlkwd{log}\hlstd{(p}\hlopt{/}\hlkwd{mean}\hlstd{(x))}\hlopt{-}\hlstd{n}\hlopt{*}\hlkwd{log}\hlstd{(}\hlkwd{gamma}\hlstd{(p))}\hlopt{+}\hlstd{(p}\hlopt{-}\hlnum{1}\hlstd{)}\hlopt{*}\hlkwd{sum}\hlstd{(}\hlkwd{log}\hlstd{(x))}\hlopt{-}\hlstd{n}\hlopt{*}\hlstd{p}
\hlkwd{plot}\hlstd{(vp,}\hlnum{.001}\hlstd{,}\hlnum{100}\hlstd{)} \hlcom{# para comprobar comportamiento}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-5-1} 
\begin{kframe}\begin{alltt}
\hlstd{sol} \hlkwb{<-} \hlkwd{optimize}\hlstd{(vp,} \hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{100}\hlstd{),} \hlkwc{maximum}\hlstd{=}\hlnum{TRUE}\hlstd{)}
\hlstd{(pMV} \hlkwb{<-} \hlstd{sol}\hlopt{$}\hlstd{maximum)}
\end{alltt}
\begin{verbatim}
## [1] 3.997042
\end{verbatim}
\begin{alltt}
\hlstd{(aMV} \hlkwb{<-} \hlstd{pMV}\hlopt{/}\hlkwd{mean}\hlstd{(x))}
\end{alltt}
\begin{verbatim}
## [1] 1.708138
\end{verbatim}
\begin{alltt}
\hlcom{## otra forma: verosimilitud perfilada en "a"}
\hlstd{vp} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{a}\hlstd{)}
    \hlstd{n}\hlopt{*}\hlstd{a}\hlopt{*}\hlkwd{mean}\hlstd{(x)}\hlopt{*}\hlkwd{log}\hlstd{(a)}\hlopt{-}\hlstd{n}\hlopt{*}\hlkwd{log}\hlstd{(}\hlkwd{gamma}\hlstd{(a}\hlopt{*}\hlkwd{mean}\hlstd{(x)))}\hlopt{+}
        \hlstd{(a}\hlopt{*}\hlkwd{mean}\hlstd{(x)}\hlopt{-}\hlnum{1}\hlstd{)}\hlopt{*}\hlkwd{sum}\hlstd{(}\hlkwd{log}\hlstd{(x))}\hlopt{-}\hlstd{n}\hlopt{*}\hlstd{a}\hlopt{*}\hlkwd{mean}\hlstd{(x)}
\hlstd{sol} \hlkwb{<-} \hlkwd{optimize}\hlstd{(vp,} \hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{100}\hlstd{),} \hlkwc{maximum}\hlstd{=}\hlnum{TRUE}\hlstd{)}
\hlstd{(aMV} \hlkwb{<-} \hlstd{sol}\hlopt{$}\hlstd{maximum)}
\end{alltt}
\begin{verbatim}
## [1] 1.708139
\end{verbatim}
\begin{alltt}
\hlstd{(pMV} \hlkwb{<-} \hlstd{aMV}\hlopt{*}\hlkwd{mean}\hlstd{(x))}
\end{alltt}
\begin{verbatim}
## [1] 3.997044
\end{verbatim}
\begin{alltt}
\hlcom{## otra forma: maximizar conjuntamente en (p,a)}
\hlkwd{optim}\hlstd{(}\hlkwd{c}\hlstd{(pMM,aMM),}
      \hlkwa{function}\hlstd{(}\hlkwc{pa}\hlstd{)}
      \hlstd{\{}
          \hlstd{p} \hlkwb{<-} \hlstd{pa[}\hlnum{1}\hlstd{]}
          \hlstd{a} \hlkwb{<-} \hlstd{pa[}\hlnum{2}\hlstd{]}
          \hlkwd{sum}\hlstd{(}\hlkwd{dgamma}\hlstd{(x,p,a,}\hlkwc{log}\hlstd{=}\hlnum{TRUE}\hlstd{))}
      \hlstd{\},}
      \hlkwc{control}\hlstd{=}\hlkwd{list}\hlstd{(}\hlkwc{fnscale}\hlstd{=}\hlopt{-}\hlnum{1}\hlstd{))}
\end{alltt}
\begin{verbatim}
## $par
## [1] 3.998465 1.709014
## 
## $value
## [1] -29.75126
## 
## $counts
## function gradient 
##       41       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
\end{verbatim}
\end{kframe}
\end{knitrout}
Se observa que las estimaciones obtenidas maximizando la verosimilitud son similares
a las del método de los momentos.
Pueden preferirse las máximo-verosímiles ya que el estimador de máxima verosimilitud
suele tener mejores propiedades que el del método de los momentos.
\item \apcua

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{## bootstrap básico paramétrico}
\hlstd{distri} \hlkwb{<-} \hlkwd{replicate}\hlstd{(}\hlnum{10000}\hlstd{,}
\hlstd{\{}
    \hlstd{mu} \hlkwb{<-} \hlkwd{mean}\hlstd{(x); s2} \hlkwb{<-} \hlkwd{var}\hlstd{(x); p} \hlkwb{<-} \hlstd{mu}\hlopt{^}\hlnum{2}\hlopt{/}\hlstd{s2; a} \hlkwb{<-} \hlstd{mu}\hlopt{/}\hlstd{s2}
    \hlkwd{mean}\hlstd{(}\hlkwd{rgamma}\hlstd{(n,p,a)}\hlopt{-}\hlkwd{mean}\hlstd{(x))}
\hlstd{\})}
\hlkwd{mean}\hlstd{(x)}\hlopt{-}\hlkwd{quantile}\hlstd{(distri,}\hlkwd{c}\hlstd{(}\hlnum{.975}\hlstd{,}\hlnum{.025}\hlstd{),}\hlkwc{names}\hlstd{=}\hlnum{FALSE}\hlstd{)}
\end{alltt}
\begin{verbatim}
## [1] 1.772318 2.858683
\end{verbatim}
\begin{alltt}
\hlcom{## bootstrap básico no paramétrico}
\hlstd{distri} \hlkwb{<-} \hlkwd{replicate}\hlstd{(}\hlnum{10000}\hlstd{,} \hlkwd{mean}\hlstd{(}\hlkwd{sample}\hlstd{(x,,}\hlnum{TRUE}\hlstd{)}\hlopt{-}\hlkwd{mean}\hlstd{(x)))}
\hlkwd{mean}\hlstd{(x)}\hlopt{-}\hlkwd{quantile}\hlstd{(distri,}\hlkwd{c}\hlstd{(}\hlnum{.975}\hlstd{,}\hlnum{.025}\hlstd{),}\hlkwc{names}\hlstd{=}\hlnum{FALSE}\hlstd{)}
\end{alltt}
\begin{verbatim}
## [1] 1.794875 2.840000
\end{verbatim}
\begin{alltt}
\hlcom{## t-bootstrap: error típico de estima mu = raíz(sigma2/n)}
\hlstd{distri} \hlkwb{<-} \hlkwd{replicate}\hlstd{(}\hlnum{10000}\hlstd{,}
                    \hlkwd{mean}\hlstd{((xb}\hlkwb{<-}\hlkwd{sample}\hlstd{(x,,}\hlnum{TRUE}\hlstd{))}\hlopt{-}\hlkwd{mean}\hlstd{(x))}\hlopt{/}\hlkwd{sqrt}\hlstd{(}\hlkwd{var}\hlstd{(xb)}\hlopt{/}\hlstd{n))}
\hlkwd{mean}\hlstd{(x)}\hlopt{-}\hlkwd{quantile}\hlstd{(distri,}\hlkwd{c}\hlstd{(}\hlnum{.975}\hlstd{,}\hlnum{.025}\hlstd{),}\hlkwc{names}\hlstd{=}\hlnum{FALSE}\hlstd{)}\hlopt{*}\hlkwd{sqrt}\hlstd{(}\hlkwd{var}\hlstd{(x)}\hlopt{/}\hlstd{n)}
\end{alltt}
\begin{verbatim}
## [1] 1.834237 3.030428
\end{verbatim}
\begin{alltt}
\hlcom{## bootstrap percentil}
\hlstd{distri} \hlkwb{<-} \hlkwd{replicate}\hlstd{(}\hlnum{10000}\hlstd{,} \hlkwd{mean}\hlstd{(}\hlkwd{sample}\hlstd{(x,,}\hlnum{TRUE}\hlstd{)))}
\hlkwd{quantile}\hlstd{(distri,}\hlkwd{c}\hlstd{(}\hlnum{.025}\hlstd{,}\hlnum{.975}\hlstd{),}\hlkwc{names}\hlstd{=}\hlnum{FALSE}\hlstd{)}
\end{alltt}
\begin{verbatim}
## [1] 1.830 2.905
\end{verbatim}
\end{kframe}
\end{knitrout}
El método paramétrico aprovecha el conocimiento de la distribución, por lo que sería en principio el mejor de los cuatro métodos.
También el $t$-\emph{bootstrap} no paramétrico, en parte,
ya que aprovecha la expresión paramétrica del error típico.  Podrían combinarse ambos:
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{## t-bootstrap paramétrico}
\hlstd{mu} \hlkwb{<-} \hlkwd{mean}\hlstd{(x); s2} \hlkwb{<-} \hlkwd{var}\hlstd{(x); p} \hlkwb{<-} \hlstd{mu}\hlopt{^}\hlnum{2}\hlopt{/}\hlstd{s2; a} \hlkwb{<-} \hlstd{mu}\hlopt{/}\hlstd{s2}
\hlstd{distri} \hlkwb{<-} \hlkwd{replicate}\hlstd{(}\hlnum{10000}\hlstd{,}
                    \hlkwd{mean}\hlstd{(((xb}\hlkwb{<-}\hlkwd{rgamma}\hlstd{(n,p,a))}\hlopt{-}\hlkwd{mean}\hlstd{(x))}\hlopt{/}\hlkwd{sqrt}\hlstd{(}\hlkwd{var}\hlstd{(xb)}\hlopt{/}\hlstd{n)))}
\hlkwd{mean}\hlstd{(x)}\hlopt{-}\hlkwd{quantile}\hlstd{(distri,}\hlkwd{c}\hlstd{(}\hlnum{.975}\hlstd{,}\hlnum{.025}\hlstd{),}\hlkwc{names}\hlstd{=}\hlnum{FALSE}\hlstd{)}\hlopt{*}\hlkwd{sqrt}\hlstd{(}\hlkwd{var}\hlstd{(x)}\hlopt{/}\hlstd{n)}
\end{alltt}
\begin{verbatim}
## [1] 1.840667 3.045840
\end{verbatim}
\end{kframe}
\end{knitrout}
\end{enumerate}
\end{enumerate}
\end{document}
