# How to calculate the Bayesian posterior analytically and by simulation?

I am working with this model:

Prior: $P(\lambda)$~ N(0, 1), only the positive part

likelihood: $P(x) = 1 - e^{-\lambda x}$ or $P(\vec{x}|\lambda)=\prod(1-e^{-\lambda x})$

Posterior: $P(\lambda|\vec{x})$

What is $P(\lambda|\vec{x})$ ?

I have something looks like this:

My understanding is that $P(x)$ is always a density function.
In here, my model is NOT a density, it is a exponential cdf. So, I am a bit confused. Is my model wrong, is this possible? Thanks!

• I always thought it is the pdf, but my colleague told me his model is $1-e^{\lambda x}$, and thus the likelihood. I think it is impossible. Could you confirm? Dec 7, 2014 at 20:36
• @Xi'an: I disagree. It makes perfect sense to use the CDF in the likelihood - just try to find the likelihood to the following question: 20 identical components were installed in an electrical substation 1000 days ago. You are asked to estimate the distribution of the time to failure of such components based on the observations of failure times in this substation. The observed values are in increasing order [in days]: 64, 396, 662, 714, 725, 848, 897, 985. 12 components are still working after 1000 days. Assume that time to failure has the Exponential distribution. Dec 8, 2014 at 10:41
• @Summit: in your censored example, the observed random variable is mixed, with pdf $$\lambda\exp\{-x\lambda\}\mathbb{I}(x<1000) + \exp(-\lambda 1000)\mathbb{I}(x=1000)$$ so even in this case one uses the pdf of the observed random variable. Dec 8, 2014 at 12:40
• @Xi'an: You are right. I was a bit sloppy in my expression. The point is: $\exp(-\lambda1000)$ is actually 1-CDF(of Exponential distribution). Imagine that the likelihood in my example is only defined with respect to the second term. However, in your answer you transformed it to the type appearing in the first term. (does this make sense to you?) Dec 8, 2014 at 12:58
• @Summit: yes indeed, the originator of the data should clarify this point. Note that a difficulty with using the version $$\prod_{i=1} [1-\exp(-\lambda x_i)]$$ is that the MLE is $\hat\lambda=+\infty$, i.e., it does not exist. Dec 8, 2014 at 14:01

A model distribution and a prior distribution are defined by either their pdf or their cdf. So your dataset is an exponential sample with parameter $\lambda$, which means that $$F_\lambda(x)=1-e^{-\lambda x}\ \text{ and }\ f_\lambda(x)=\lambda e^{-x\lambda}$$ are the cdf and pdf for this model. Therefore the likelihood function is defined as $$L(\lambda|x_1,\ldots,x_n) = \prod_{i=1}^n f_\lambda(x_i) = \lambda^n \exp\left\{-\lambda\sum_{i=1}^n x_i\right\}\,.$$ If your prior on $\lambda$ is a truncated normal $\text{N}(0,1)$, then its pdf is given by $$\pi(\lambda) = \sqrt{\dfrac{2}{\pi}}\ \exp\{-\lambda^2/2\}\mathbb{I}_{\lambda>0}\,.$$ (The cdf is irrelevant for the computation of the posterior.)

Now the posterior distribution is defined via Bayes' formula by $$\pi(\lambda|x_1,\ldots,x_n) \propto \pi(\lambda) L\lambda|x_1,\ldots,x_n)\,.$$ Hence, \begin{align*} \pi(\lambda|x_1,\ldots,x_n) &\propto \lambda^n \exp\{-\lambda^2/2\} \exp\left\{-\lambda\sum_{i=1}^n x_i\right\}\mathbb{I}_{\lambda>0}\\ &\propto \lambda^n \exp\left\{-\frac{1}{2}\left(\lambda+\sum_{i=1}^n x_i\right)^2\right\}\mathbb{I}_{\lambda>0} \end{align*} which is the density of a non-standard distribution.

• (+1) I wonder why he would use this prior when a conjugate one is available. He can even use a mixture of gamma's as his prior, which is quite flexible.
– Zen
Dec 7, 2014 at 20:53
• Thanks for the detailed explanation. Is it common practice take log likelihood? I can never accurately decide. Your work looks great without it. Dec 7, 2014 at 21:34
• The last sentence "which is the density of a non-standard distribution." My colleague wants me to give the distribution of $\lambda$, he wants the form $1-e^{-\lambda x}$, is that possible? I don't know how to get it. Dec 7, 2014 at 22:56
• Summit wonders if your data is made from exponential $x_i\sim\mathcal{E}(\lambda)$ observations or from censored exponential durations $x_i$, ie you only know that the exponential durations are at most $x_i$, which happens with probability $1-e^{-\lambda x_i}$. Dec 8, 2014 at 12:53
• @Xi'an I discussed with my colleague today. The model is $y = 1- e^{-\lambda{x}}+\epsilon$ where, $\epsilon \sim N(0, 1)$. So, the likelihood model the product of $y_i \sim N(1-e^{\lambda{x_i}}, 1)$. Dec 8, 2014 at 22:00

Suppose that $X_1,\dots,X_n$, given $\Lambda=\lambda$, are conditionally independent and identically distributed, such that $X_1\sim\mathrm{Exp}(\lambda)$, and you believe a priori that $\Lambda\sim\mathrm{Ga}(\alpha_0,\beta_0)$, in which $\alpha_0,\beta_0>0$ are specified real numbers.

After you get a sample $x=(x_1,\dots,x_n)$, the likelihood is the sampling density seem as a function of $\lambda$. $$L_x(\lambda) = f_{X_1,\dots,X_n\mid\Lambda}(x_1,\dots,x_n\mid\lambda) = \prod_{i=1}^n f_{X_i\mid\Lambda}(x_i\mid\lambda) = \lambda^n e^{-\lambda\sum_{i=1}^n x_i} \,I_{(0,\infty)}(\lambda)\, .$$

Bayes' Theorem gives $$f_{\Lambda\mid X_1,\dots,X_n}(\lambda\mid x_1,\dots,x_n) \propto L_x(\lambda)\,f_\Lambda(\lambda) \propto \lambda^{\alpha_0+n-1} e^{-(\beta_0+\sum_{i=1}^n x_i)\lambda}\,I_{(0,\infty)}(\lambda) \, ,$$ yielding that a posteriori $$\Lambda\mid X_1=x_1,\dots,X_n=x_n\sim\mathrm{Ga}\!\left(\alpha_0+n,\beta_0+\sum_{i=1}^n x_i\right) \, .$$ The Bayes estimate with quadratic loss is $$\mathrm{E}\left[\Lambda\mid X_1=x_1,\dots,X_n=x_n\right] = \frac{\alpha_0+n}{\beta_0+\sum_{i=1}^n x_i} \, .$$ A credible interval can be obtained by Monte Carlo.

a0 <- 6
b0 <- 20

n <- 45
x <- rexp(n, rate = 5)

N <- 10^6

level <- 0.95

post <- rgamma(N, shape = a0 + n, rate = b0 + sum(x))

cat("Bayes estimate:", (a0+n)/(b0+sum(x)), "\n")
cat(sprintf("%.2f", level*100), "% Credible Interval: [ ",
quantile(post,(1-level)/2), " ; ",
quantile(post,(1+level)/2), " ]\n", sep = "")

• Then I'm afraid your model makes no sense. For each fixed $\lambda$, the likelihood must be a density/pmf for the $x_i$'s. Try to check this in your case.
– Zen
Dec 7, 2014 at 20:22
• You have to apply the correct version of Bayes' formula: $$\pi(\lambda|x) \propto \pi(\lambda) p(x|\lambda)$$ operates on densities (pdfs), not cdfs. And the likelihood of a sample is the product of the individual densities, not of the individual cdfs. Dec 7, 2014 at 20:30
• @user13985: you're rushing and trying to do inference without understanding the concepts. This book is an excellent option: tinyurl.com/kqsnycq.
– Zen
Dec 7, 2014 at 20:36
• Ask him to consider $\lambda=1$ and sum over $x$ (just one observation will be enough).
– Zen
Dec 7, 2014 at 20:39
• Thanks for buying the book, you should re-check Chapter 1 for the proper definitions of the prior and likelihood. Dec 7, 2014 at 20:51