Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

! Package amsmath Error: \hat allowed only in math mode. Error: Failed to compile MLE.tex. See MLE.log for more info. Execution halted #4089

Closed
zjplab opened this issue Dec 26, 2018 · 8 comments
Assignees
Labels

Comments

@zjplab
Copy link

zjplab commented Dec 26, 2018

System details

[1] "1.1.456"

$R
[1] "C:\PROGRA1\MICROS2\ROPEN1\R-351.1\bin\x64\R.exe"

$pdflatex
[1] "C:\Users\zjplab\AppData\Roaming\TinyTeX\bin\win32\pdflatex.exe"

$bibtex
[1] "C:\Users\zjplab\AppData\Roaming\TinyTeX\bin\win32\bibtex.exe"

$gcc
[1] "C:\Users\zjplab\MINICO~1\Library\MINGW-~1\bin\gcc.exe"

$git
[1] "C:\PROGRA~1\Git\cmd\git.exe"

$svn
[1] ""

R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] RevoUtils_11.0.1 RevoUtilsMath_11.0.0

loaded via a namespace (and not attached):
[1] compiler_3.5.1 backports_1.1.2 magrittr_1.5 rprojroot_1.3-2 htmltools_0.3.6
[6] tools_3.5.1 yaml_2.2.0 Rcpp_0.12.18 stringi_1.1.7 rmarkdown_1.10
[11] knitr_1.20 stringr_1.3.1 digest_0.6.15 evaluate_0.11

SysInfo:
sysname release version nodename machine
"Windows" ">= 8 x64" "build 9200" "WINDOWS-IO3NS2G" "x86-64"
login user effective_user
"zjplab" "zjplab" "zjplab"

R Version:
_
platform x86_64-w64-mingw32
arch x86_64
os mingw32
system x86_64, mingw32
status
major 3
minor 5.1
year 2018
month 07
day 02
svn rev 74947
language R
version.string R version 3.5.1 (2018-07-02)
nickname Feather Spray
$ALLUSERSPROFILE
[1] "C:\ProgramData"

Steps to reproduce the problem

Describe the problem in detail

Describe the behavior you expected

You could try knit the following RNotebook Code to pdf to reproduce this error?

---
title: "Maximum Likehood Estimator"
output:
  pdf_document: default
  html_notebook: default
  word_document: default
  html_document:
    df_print: paged
---
# Introduction

maximum likelihood estimation is one common method for estimating paremater in a parametric model, just like method of moments. We assume$X_1, X_2..., X_n$ be IID with *PDF* $f(x;\theta)$, define the likelihood function as 
$$
\mathcal{L}(\theta)=\Pi_{i=1}^nf(X_i;\theta)
$$
And the log likelihood function is defined by $\mathcal{l}_n(\theta)=log \mathcal{L}_n(\theta)=\sum_{i=1}^n log f(X_i; \theta)$. 


The *maximum likelihood estimator* MLE, denoted by $\hat{\theta}_n$, is the value of $\theta$ that maximizes $\mathcal{L}_n(\theta)$




# Derivation 


### Point mass at $a$
Notice the probability mass function for point mass distribution is $f(x)=1$ in $x=a$ and 0 elsewhere. In fact, $f$ have only two values:1 and 0, so does the $\mathcal{L}(\theta)$ which is the product of many $f$s. We choose the MLE for point mass distribution to be mode among $X_n$. 

### Bernoulli
The probability function is $f(x;p)=p^x(1-p)^{1-x}$ for $x=0,1$ , the unknown parameter is $p$.
$$
\mathcal{L}_n(p)=\prod_{i=1}^nf(X_i; p)=\prod_{i=1}^n p^{X_i}(1-p)^{1-X_i}=p^S(1-p)^{n-S}
$$
where $S=\sum_i X_i$. Hence,
$$
\mathcal{l}_n(p)=S log p+ (n-S)log(1-p)
$$
Take the derivation and set it equal to 0 gave us $\hat{p}_n=\frac{S}{n}$

### Binomial(N,p)
$$
\mathcal{L}(p) = \prod_i^n(f(y_i)) =  \prod_i^n \left[ {{N}\choose{y_i}}p^{y_i} (1-p)^{N- yi} \right] = \left[ \prod_i^n {{N}\choose{y_i}} \right] p^{\sum_1^n y_i} (1-p)^{nN - \sum_1^n{y_i}} \,.$$


$$\hat{p} = \arg\max_p \left[ \left[ \prod_i^n {{N}\choose{y_i}} \right] p^{\sum_1^n y_i} (1-p)^{nN - \sum_1^n{y_i}}\right]$$

$$
\begin{array}{rcl} \hat{p} & = & \displaystyle\arg\max_p\ \  \ln\left[ \left[ \prod_i^n {{N}\choose{y_i}} \right] p^{\sum_1^n y_i} (1-p)^{nN - \sum_1^n{y_i}}\right] \\ & = & \displaystyle \arg\limits\max_p  \sum_{i=1}^n\left[ \ln {{N}\choose{y_i}} + \left( \sum_{i=1}^n y_i \right)\ln(p) + \left( nN - \sum_{i=1}^n y_i \right)\ln(1-p) \right] \end{array}
$$

$$
0 + \frac{\left( \sum_{i=1}^n y_i \right)}{\hat{p}} - \frac{\left(n N - \sum_{i=1}^n y_i \right)}{1-\hat{p}} = 0 \,.
$$
$$\hat{p} = \frac{\sum_{i=1}^n y_i}{nN} = \frac{\bar{y}}{N} $$

### Geometric(p)
$P(X=k)=p(1-p)^k$ where $k\ge1$
$$
\mathcal{L}\left(p \right)={\left(1-p \right)}^{{x}_{1}-1}p {\left(1-p \right)}^{{x}_{2}-1}p...{\left(1-p \right)}^{{x}_{n}-1}p ={p}^{n}{\left(1-p \right)}^{\sum_{1}^{n}{x}_{i}-n}
$$
$$ln\mathcal{L}(p)=\sum_{i=1}^{n}{x}_{i}ln(p)+\left(n-\sum_{i=1}^{n}{x}_{i} \right)ln\left(1-p \right)$$

$$
\frac{dln\mathcal{L}(p)}{dp}=\frac{1}{p}\sum_{i=1}^{n}{x}_{i}+\frac{1}{1-p}\left(n-\sum_{i=1}^{n}{x}_{i} \right)=0
$$
$$\hat{p}=\frac{n}{\left(\sum_{1}^{n}{x}_{i} \right)}$$

### Poisson Distribution($\lambda$)
$$
\mathcal{L}(\lambda)=\prod_{i=1}^{n}\frac{{\lambda}^{{x}_{i}}{e}^{-\lambda}}{{x}_{i}!} = {e}^{-n\lambda} \frac{{\lambda}^{\sum_{1}^{n}{x}_{i}}}{\prod_{i=1}^{n}{x}_{i}!}
$$
$$
lnL(\lambda)=-n\lambda+\sum_{1}^{n}{x}_{i}ln(\lambda)-ln\left(\prod_{i=1}^{n}{x}_{i}!\right)
$$
$$\frac{dlnL(\lambda)}{dp}=-n+   \sum_{1}^{n} \frac{ x_i } {\lambda} $$

$$\hat{\lambda}=\frac{\sum_{i=1}^{n}{x}_{i}}{n}$$

### Uniform(a,b)
$$
\mathcal{L}_n(a, b)=\frac{1}{(b-a)^n}
$$
if all $a \le X_1, X_2, ..., X_n \le b$ 
else $\mathcal{L}=0$
To maximize, we need to minize $b-a$, the boundary condition is to choose $\hat{a}=min_i(X_i)$ and $\hat{b}=max_i(X_i)$

### Normal($\mu$, $\sigma^2$)
$$
\mathcal{L}_n(\mu, \sigma) =\prod_i \frac{1}{\sigma} exp\{- \frac{1}{2 \sigma^2} (X_i-\mu)^2 \} \\
                          =\sigma^{-n}exp\{ -\frac{1}{2\sigma^2}\sum_i(X_i-\mu)^2 \}\\
                          =\sigma^{-n}exp\{ -\frac{nS^2}{2\sigma^2}\} exp\{- \frac{n(\bar{X} -\mu)^2}{2\sigma^2} \}
$$
where $\bar{X}=n^{-1}\sum_i X_i$, $S^2=n^{-1}\sum_i (X_i - \bar{X})^2$ 

$$
\mathcal{l}(\mu, \sigma)=-n log\sigma - \frac{nS^2}{2\sigma^2}- \frac{n(\bar{X}-\mu )^2}{2\sigma^2}
$$
Solving eqations $\frac{ \partial \mathcal{l}(\mu, \sigma) }{\partial \mu} =0$ and 
$\frac{ \partial \mathcal{l}(\mu, \sigma)} {\partial \mathcal{l} \sigma}$
we conclude that $\hat{\mu}=\bar{X}$ and $\hat{\sigma}=S$

### Exponential($\beta$)
$$
\mathcal{L}(\beta)=\mathcal{L}\left(\beta;{X}_{1},{X}_{2}...{X}_{n} \right)=\left(\frac{1}{\beta}{e}^{\frac{{-X}_{1}}{\beta}}\right)\left(\frac{1}{\theta}{e}^{\frac{{-X}_{2}}{\beta}}\right)...\left(\frac{1}{\theta}{e}^{\frac{{-X}_{n}}{\beta}} \right)=\frac{1}{{\beta}^{n}}exp\left(\frac{-\sum_{1}^{n}{X}_{i}}{\beta} \right)
$$
$$
ln \mathcal{L}\left(\beta\right)=-n ln\left(\beta\right) -\frac{1}{\beta}\sum_{1}^{n}{X}_{i}, 0<\beta<\infty$$
$$
\frac{d\left[ln \mathcal{L}\left(\theta\right) \right]}{d\beta}=\frac{-n }{\beta} +\frac{1}{{\beta}^{2}}\sum_{1}^{n}{X}_{i}=0
$$

$$
\hat{\beta}=\frac{\sum_{1}^{n}{X}_{i}}{n}
$$

### Gamma($\alpha, \beta$)
$$
\mathcal{L}(\alpha, \beta)=(\frac{1}{\beta^\alpha \Gamma(\alpha)})^n exp^{-\frac{\sum_iX_i}{\beta}} \prod_{i=1}^n X_i^{\alpha-1}
$$



$$
\mathcal{l}_n(\alpha, \beta)=
-n \alpha log \beta -n log\Gamma(\alpha) -\frac{\sum_i X_i}{\beta}
+\sum_{i=1}^n (\alpha -1) log X_i
$$
Solve $\frac{\partial \mathcal{l}_n(\alpha, \beta)}{\partial \beta}=0$ and $\frac{\partial l_n(\alpha, \beta)}{ \partial \alpha}=0$  yield:
$$
\hat{\beta}=\frac{\sum_i X_i}{n\alpha}
$$

and $$
\sum_{i=1}^n log X_i = n og \frac{\sum X_i}{n \alpha} + n \psi(\alpha) 
$$

where $$\psi(\alpha)=\frac{d log \Gamma(\alpha) }{d \alpha}$$, is called Digamma Function
However, no closed form exists for $\hat{\alpha}$, only numerical solution exist. 

we calculate the numerical estimation using Newton-Raphson [Refernce](http://bioops.info/2015/01/gamma-mme-mle/)
$$
l_n'(\alpha)=n log(\frac{\alpha}{\bar{X}})-n\frac{\Gamma'(\alpha)}{\Gamma(\alpha)}+\sum_{i=1}^nlog X_i
$$
$$
l_n''(\alpha)=\frac{n}{\alpha}-n(\frac{\Gamma'(\alpha)}{\Gamma(\alpha)})'
$$
$$
\hat{\alpha}^{(k)}=\hat{\alpha}^{(k-1)}- \frac{l_n'(\hat{\alpha} ^{(k-1)})} {l_n''(\hat{\alpha} ^{(k-1)})}
\\
\hat{\beta}=\frac{\bar{X}}{\hat{\alpha}}
$$
### Beta($\alpha, \beta$)
$$
\mathcal{L}_n(\alpha, \beta)=(\frac{\Gamma(\alpha +\beta)}{\Gamma(\alpha)\Gamma(\beta)})^n \prod_{i=1}^n(X_i)^{\alpha-1}
\prod_{i=1}^n(1-X_i)^{\beta-1}
$$
$$
\mathcal{l}_n(\alpha, \beta)=nlog( \Gamma (\alpha +\beta) ) -n log( \Gamma(\beta)) +(\alpha-1)\sum_{i=1}^n log(X_i)+(\beta-1) \sum_{i=1}^n log(1-X_i)
$$
solve $\frac{\partial \mathcal{l}_n(\alpha, \beta)}{ \partial \alpha}$ and 
$\frac{ \partial \mathcal{l}_n(\alpha, \beta)}{\partial \beta}$
 yield no analytical solution. 
 However, numerical solution via Newton-Raphson mehod does exist.  [For reference, check page 27 on this paper](https://scholarsarchive.byu.edu/cgi/viewcontent.cgi?article=2613&context=etd)

We set $\hat{\theta}=(\hat{\alpha}, \hat{ \beta})$ iteratively:
$$
\hat{\theta}_{i+1}=\hat{\theta}_{i}- G^{-1}g\\
g=[g1,g2]\\
g1=\psi(\hat{\alpha})-\psi(\hat{\alpha}+\hat{\beta})-\frac{1}{n}\sum_{i=1}^n log(X_i)\\
g2=\psi(\hat{\beta})-\psi(\hat{\alpha}+\hat{\beta})-\frac{1}{n}\sum_{i=1}^n log(1-X_i)\\
G=\begin{bmatrix}
\frac{dg1}{d\alpha} & \frac{dg2}{d\beta} \\
\frac{dg2}{d\alpha} & \frac{dg2}{d\beta}

\end{bmatrix}\\
\frac{dg1}{d\alpha}=\psi'(\alpha) - \psi'(\alpha +\beta) \\
\frac{dg1}{d\beta}=\frac{dg2}{d\beta}=-\psi'(\alpha + \beta)\\
\frac{dg2}{d\beta}=\psi'(\beta)-\psi'(\alpha+\beta)
$$
 
### Student t distribution
$$
\mathcal{L}_n(v)=(\frac{\Gamma( \frac{\nu+1}{2} )}{\Gamma{\frac{ \nu}{ 2}} })^n \frac{1}
{
\prod_{i=1}^n(1+\frac{X_i^2}{\nu})^{ \frac{\nu+1} {2} }
} 
$$
$$
l_n(\nu)=n log(\Gamma(\frac{\nu+1}{2}))-n log( \Gamma( \frac{\nu}{2} ) )-\sum_{i=1}^n \frac{\nu+1}{2}log(1+ \frac{X_i^2}{\nu} )
$$
$$\frac{\partial l_n(\nu)}{\partial\nu}= \frac{\nu+1}{2}\sum_{i=1}^n \frac{X_i^2}{\nu^2}-\frac{1}{2}\sum_{i=1}^n(\frac{X_i^2}{\nu}+1)+\frac{n \psi( \frac{\nu+1}{2})}{2}- \frac{n \psi(\frac{\nu}{2})}{2}=0
$$
We cannot get closed-form solution. 

### $\mathcal{\chi}_p^2$ Chi-Squared Distriubtion
$$
\mathcal{L}_n(p)=\frac{1}{(\Gamma(p/2) 2^{p/2})^n}e^{\frac{-\sum_{i=1}^nX_i}{2}} \prod_{i=1}^n X_i^{p/2-1}
$$

$$
l_n(p)=-n log(\Gamma(p/2))- \frac{pn}{2} log2-\sum_{i=1}^n\frac{X_i}{2}+(\frac{p}{2}-1)\sum_{i=1}^n log X_i
$$

$$
\frac{\partial l_n(p)}{\partial p}=
- \frac{n}{2}\psi(p/2)-\frac{n}{2} log2+\frac{1}{2}\sum_{i=1}^n log X_i=0
$$

$$
\hat{p}=2 \psi^{-1}(\frac{\sum_{i=1}^n log X_i}{n} - log 2)
$$
But [in fact](http://mathworld.wolfram.com/DigammaFunction.html), digamma function is not monotone. Therefore it doesn't have an inverse function at all. So the formula above is not right! In other words, the student t distributions doesn't have a closed-form expression for MLE. 
@ronblum
Copy link
Contributor

ronblum commented Dec 26, 2018

@zjplab Line 175 is blank. If you remove it, the file should knit properly.

Change

\frac{dg2}{d\alpha} & \frac{dg2}{d\beta}

\end{bmatrix}\\

to

\frac{dg2}{d\alpha} & \frac{dg2}{d\beta}
\end{bmatrix}\\

@ronblum ronblum closed this as completed Dec 26, 2018
@ronblum ronblum added info needed Additional information requested—reprex, steps, open question, etc. and removed info needed Additional information requested—reprex, steps, open question, etc. labels Dec 26, 2018
@zjplab
Copy link
Author

zjplab commented Dec 28, 2018

@ronblum Thanks. But just out of curiosity, why does a blank lead to such error?

@ronblum
Copy link
Contributor

ronblum commented Dec 28, 2018

@zjplab Empty lines are not allowed in Latex math mode, according to this LaTeX wikibook and a thread on Stack Exchange.

@brainfo
Copy link

brainfo commented Mar 18, 2020

Hi, I've encountered the same error info:

! Package amsmath Error: \hat allowed only in math mode.

错误: LaTeX failed to compile Bayes-Assignment-2.tex. See https://yihui.org/tinytex/r/#debugging for debugging tips. See Bayes-Assignment-2.log for more info.
停止执行

and I've removed all empty lines:

---
title: "Bayes Assignment Chapter3"
author: "Hong JIANG"
date: "2020年3月11日"
#output: pdf_document
documentclass: ctexart
output:
  rticles::ctex:
    fig_caption: yes
    number_sections: no
    toc: yes
classoption: "hyperref"
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
Sys.setlocale(category="LC_ALL",locale="en_US.UTF-8") 

3.11

证:

(1)

左式:
令$s^2=\frac{1}{n-1}\sum_{i=1}^{x}(X_i-\bar{X})^2$,$T=(\bar{X},s^2)$,$\varphi=(\theta,\sigma^2)$,
则有$\bar{X}|\varphi\sim N(\theta,\frac{\sigma^2}{n})$,$\frac{(n-1)S^2}{\sigma^2}\sim \chi^2_{n-1}$,且$\bar{X}$和$S^2$独立。
由因子分解定理,可知T是$\varphi$的充分统计量,其联合密度为
$$ f(t|\varphi)=c\cdot \sigma^{-n}(s^2)^{\frac{n-1}{2}-1}e^{-\frac{1}{2\sigma^2}[(n-1)s^2+n(\bar{x}-\theta)^2]}$$
其中$c=\frac{\sqrt{n}(n-1)^{\frac{n-1}{2}}}{\sqrt{\pi}\cdot 2^{\frac{n}{2}} \Gamma(\frac{n-1}{2})}$,
因而$(\theta,\sigma^2)$的似然函数
$$l(\theta,\sigma^2|t)\propto (\sigma^2)^(-\frac{n}{2})e^{-\frac{1}{2\sigma^2}}[(n-1)s^2+n(\bar{x}-\theta)^2]$$
已知$(\theta,\sigma^2)$的联合先验分布为正态-逆伽马分布,则其联合密度为:
$$
\begin{aligned}
\pi(\theta,\sigma^2)& =\pi_1(\theta|\sigma^2)\pi_2(\sigma^2) \
& \sim (\sigma^2)^(\frac{2\alpha+1}{2}+1)e^{-\frac{1}{2\sigma^2}[\frac{(\theta-\mu)^2}{\tau}+2\beta^{-1}]}
\end{aligned}
$$
则有
$$
\begin{aligned}
\pi(\theta,\sigma^2|t)&\sim l(\theta,\sigma^2|t)\pi(\theta,\sigma^2) \
&\sim (\alpha^2)^{-\frac{n+2\alpha+1}{2}+1}e^{-\frac{1}{2\alpha^2}[(n-1)s^2+n(\bar(X)-\theta)^2+\frac{(\theta-\mu)^2}{\tau}+2\beta^{-1}]} \
&\sim \sigma^2[-(\frac{n+\alpha+1}{2}+1)]e^{-\frac{1}{2\alpha^2}[(n-1)s^2+2\beta^-1+H]}
\end{aligned}
$$
其中$$
\begin{aligned}
H&=n(\bar{x}-\theta)^2+\frac{(\theta-\mu)^2}{\tau} \
&=(n+\frac{1}{\tau})(\theta-\frac{n\bar(x)+\frac{\mu}{\tau}}{n+\frac{1}{\tau}})^2+\frac{n}{\tau n+1}(\bar{x}-\mu)^2
\end{aligned}
$$.
因此,
$$ \pi(\theta,\sigma^2|t) \sim (\sigma^2)^{-[\frac{\nu_n+1}{2}+1]}e^{-\frac{1}{2\sigma_n^2}[\nu_n\sigma_n^2+\frac{(\theta-\mu(\bar{x}))^2}{\tau_n}]}$$,
此处$\nu_n=n+2\alpha$,$\tau_n=n+\frac{1}{\tau}$,$\mu(\bar{x})=\frac{\tau n \bar{x}+\mu}{n\tau+1}$,$\nu_n\sigma_n^2=(n-1)s^2+2\beta^{-1}+\frac{n}{\tau n+1}(\bar{x}-\mu)^2$.
得到后验密度为
$$\pi(\theta,\sigma^2|t)=C (\sigma^2)^{-[\frac{\nu_n+1}{2}+1]}e^{-\frac{1}{2\sigma_n^2}[\nu_n\sigma_n^2+\frac{(\theta-\mu(\bar{x}))^2}{\tau_n}]}$$,
其中$C=\sqrt{\frac{\tau_n}{2\pi}}\cdot \frac{(\nu_n\sigma_n^2)^{\frac{\nu_n}{2}}}{2^{\nu_n}{2}\Gamma(\frac{\nu_n}{2})}$.
右边:
$$
\begin{aligned}
\pi_1(\theta|\sigma^2,x)&=\frac{1}{\eta \sqrt{2\pi}}e^{-\frac{(x-\mu(x))^2}{2\sigma^2}} \
& \sim \sqrt(\frac{n+\tau^{-1}}{\sigma^2})e^{-(n+\tau^{-1})\frac{(x-\frac{\mu+n\tau \bar{x}}{n\tau+1})^2}{2\sigma^2}}
\end{aligned}
$$
$$\pi_2(\sigma^2,x)=\frac{e^{-\frac{\bar{\beta}}{x}}(\frac{\bar{\beta}}{x})^{\alpha+\frac{n}{2}}}{x\Gamma(\alpha+\frac{n}{2})}$$
$\therefore$ $\pi_1(\theta|\sigma^2,x)\cdot \pi_2(\sigma^2,x)=\frac{1}{\eta \sqrt{2\pi}}\cdot \frac{1}{x\Gamma(\alpha+\frac{n}{2})}\cdot (\frac{\bar{\beta}}{x})^(\alpha+\frac{n}{2})$.
得证。

(2)

对联合后验密度关于$\theta$积分,有:
$$
\begin{aligned}
\pi(\sigma^2|t)&=\int \pi(\theta,\sigma^2|t) d\theta \
&= \frac{(\nu_n \sigma_n^2)^{\frac{\nu_n}{2}}}{2^{\frac{\nu_n}{2}}\Gamma(\frac{\nu_n}{2})} \cdot (\sigma^2)^{-(\frac{\nu_n}{2}+1)} \cdot e^{-\frac{\nu_n\sigma_n^2}{2\sigma^2}}
\end{aligned}
$$
由于$\nu_n=n+2\alpha$,$\tau_n=n+frac{1}{\tau}$,
该边缘密度函数为$\Gamma^{-1}(\nu_n/2,\nu_n\tau_n^2/2)$,的密度函数即为$\Gamma^{-1}(\alpha+n/2,\widetilde{\beta})$的密度函数。

(3)

对联合后验密度关于$\sigma^2$积分,有:
$$
\begin{aligned}
\pi(\theta|t)&=\int \pi(\theta,\sigma^2|t) d\sigma^2 \
&= C 2^{\frac{\nu_n+1}{2}}\Gamma(\frac{\nu_n+1}{2})[\nu_n\sigma_n^2+\tau_n(\theta-\mu(\bar{x})^2)]^{-\frac{\nu_n+1}{2}} \
&= \frac{\Gamma(\frac{\nu_n+1}{2})}{\Gamma(\frac{\nu_n}{2})\sqrt(\pi \nu_n)} \cdot \frac{\sqrt(\tau_n)}{\sigma_n}[1+\frac{1}{\nu_n}(\frac{\theta-\mu(\bar{x})}{\sigma_n/\sqrt{\tau_n}})^2]^{-\frac{\nu_n+1}{2}}
\end{aligned}
$$
这是一元t分布$t(\nu_n,\mu(\bar{x}),\frac{\sigma_n^2}{\tau_n})$的密度函数。替换$\nu_n$和$\tau_n$得证。

3.15

(1)

解:只对X做一次观察时:
$$
\begin{aligned}
h(x,\theta)&=f(x|\theta)\pi(\theta) \
&=\theta (1-\theta)^(x-1), 0\leq \theta \leq 1
\end{aligned}
$$
$$
\begin{aligned}
m(x)&=m(x|\pi)=\int_\Theta f(x|\theta)\pi(\theta)d\theta \
&= \int_0^1 \theta(1-\theta)^{x-1}d\theta \
&= -\int_0^1 (1-t)t^{x-1}dt \
&= \left[\frac{1}{x+1}t^{x+1}\right]_0^1 - \left[\frac{1}{x}t^x\right]_0^1 \
&= -\frac{1}{x(x+1)}
\end{aligned}
$$

则$\pi(\theta|x)=\frac{h(x,\theta)}{m(x)}=-x(x+1)\theta(1-\theta)^{x-1}$.
有$$
\begin{aligned}
\hat{\theta_E}=E(\theta|x)&=-x(x+1)\int_0^1 \theta^2(1-\theta)^{x-1}d\theta \
&= \int_0^1 (1-t)^2t^(x-1)dt \
&= 2\frac{x+1}{x+2}
\end{aligned}
$$
当$X=3$,有$\hat{\theta_E}=\frac{8}{5}$.

(3)

$f(X_1=3,X_2=2,X_3=5,\theta)=\theta^3(1-\theta)^7$
$$
m(x|\pi)=\int_0^1 \theta^3(1-\theta)^7
$$

$\pi(\theta|x)=\frac{f(X_1=3,X_2=2,X_3=5,\theta)}{m(x|\pi)}$
则$$\hat{\theta|X_1=3,X_2=2,X_3=5}=\int_0^1 \theta \pi(\theta|x) d\theta $$

library(MASS)
integrand <- function(theta){
  theta^3*(1-theta)^7
}
integrand2 <- function(theta){
  theta^4*(1-theta)^7
}
thetaHat <- fractions(integrate(integrand2,lower = 0,upper = 1)$value/integrate(integrand,lower = 0,upper = 1)$value)

解得$\hat{\theta_E}=1/3}

3.28

(1)

$\theta$的后验分布$\pi(\theta|x)$为$(\mu_n(\bar{X}),\eta_n^2)$,其中
$$
\begin{aligned}
\mu_n(\bar{X})&=\frac{2.19}{\frac{1}{n}+2.19}\bar{X}=6\cdot \frac{2.19}{3.19} \
\eta_n^2 &= \frac{2.19}{1+2.19n} = \frac{2.19}{3.19}
\end{aligned}
$$
$\because$ 后验分布单峰对称,$\therefore$ HPD可信区间为等尾区间。
令$$P(\mu_n(\bar{X}))-\eta_n u_{\frac{\alpha}{2}} \leq \theta \leq \mu_n(\bar{X}))+\eta_n u_{\frac{\alpha}{2}}|\bar{X})=1-\alpha$$,
其中$u_{\alpha/2}$为$N(0,1)$的上策$\alpha/2$分位数,故$\theta$的$1-\alpha$的HPD可信区间为:

qnorm(c(0.05,0.95),6*2.19/3.19,sqrt(2.19/3.19))

$$\left[2.756254,5.481991\right]$$

(2)

样本$X=x_1=6$的联合分布为:
$$f(x|\theta)=\frac{1}{2\pi}e^{-\frac{(\theta-x_1)}{2}}\cdot \frac{1}{\pi(1+\theta^2)}$$
则$\theta$的后验密度为:
$$
\begin{aligned}
\pi(\theta|x_1)&=\frac{f(x_1|\theta)}{\frac{1}{\pi\sqrt{2\pi}} \int e^{-\frac{(\theta-x_1)}{2}}\cdot \frac{1}{1+\theta^2} d\theta} \
&\propto f(x|\theta)
\end{aligned}
$$
计算95%HPD可信区间:

library(GoFKernel)
library(TeachingDemos)
## posterior
dtest =function(theta) {dnorm(theta,6,1)*dcauchy(theta,0,1)}
P.density <- function(theta){
  dnorm(theta,6,1)*dcauchy(theta,0,1)/integrate(f=dtest,
                                                lower = -Inf,
                                                upper = Inf)$value
}
## CDF
CDF <- function(x){
  integrate(f=P.density,
            lower=-Inf,
            upper=x)$value
}
## calculate the inverse cdf of the posterior
## to inverse sample
icdf <- inverse(CDF)
hpd(icdf)

可信区间为$\left[ 3.941,20.76368\right]$.

3.33

有$\pi(\theta|x)\sim N(\mu_3(\bar{x}),\eta_3^2)=N(3,(1/2)^2)$

a0 <- pnorm(0.1,3,0.5)
a1 <- 1-a0

$\pi(\theta) \sim N(3,1)$, $\therefore$

pi0 <- pnorm(0.1,3,1)
pi1 <- 1-pi0
B <- a0*pi1/(a1*pi0)
S <- B>1 #FALSE

可见支持$H_0$的贝叶斯因子不高,故否定$H_0$,接受$H_1$.


Do you know how to fix it?

@EpicScizor
Copy link

@brainfo To you and someone else who later discovers via google (like I did):

It's not just newlines, but also whitespace after the dollar signs on the same line that causes this.

In @brainfo's case it's $$ f(t|\varphi)and $$ \pi which causes the error. Should be $$f(t|\varphi)and $$\pi, respectively (no whitespace after the dollar sign - only characters or newline then character)

@daniel7009
Copy link

Hi did the error leave after editing your preamble in rmarkdown ! Package amsmath Error: \hat allowed only in math mode. getting the same

@proudjiao
Copy link

proudjiao commented Jan 14, 2022

@brainfo To you and someone else who later discovers via google (like I did):

It's not just newlines, but also whitespace after the dollar signs on the same line that causes this.

In @brainfo's case it's $$ f(t|\varphi)and $$ \pi which causes the error. Should be $$f(t|\varphi)and $$\pi, respectively (no whitespace after the dollar sign - only characters or newline then character)

Thank you so much for the help. I had the same problem but was able to fix them after I delete all the white space after the single and double dollar sign!

@brainfo
Copy link

brainfo commented Jan 14, 2022 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

6 participants