# Estimating mean of Normal with unknown variance and then predict the future observation

I am trying to estimate population mean of 9 observations when the variance is unknown. I marginalized the posterior and understand that the t- distribution would give me the distribution of population mean. I am stuck at this point. Normally, If I had to estimate some thing I would generate 1000 or more random samples of the given distribution and then generate point or interval estimates for it's values. But the T-distribution has confused me. Matlab's tpdf generates only 8 samples, but when I sum them up they do not add up to 1 which looks weird so is it generating actual values? If these are actual values then where is the distribution? how do I estimate mean from it (Substitute these values in the standardization formula to find values of mean?).

PS: I have been studying stats recently and though I understand the mathematical part of it, I feel miserable when doing simulation in matlab. So I would appreciate any pointers twards learning the computational side of it.

EDIT: I understand the mathematical or derivation part of it. It is the computational simulation that confuses me. I use tpdf for using t distribution but it needs data and degree of freedom. and then how do I go about finding the point estimate of mean in matlab. Aso tpdf needs to be translated towards my data values.

• slef-study?! Please add the tag. And explain more clearly why the $t$ distribution confuses you. It has a mean that you can use as an estimator and you can produce a credible interval by using the t pdf. Apr 3, 2015 at 7:18
• Why should they add to 1? Apr 3, 2015 at 8:33
• Shouldn't a probability distribution integrate to 1?
– Sudh
Apr 3, 2015 at 15:41
• Values sampled from a distribution are not (usually) probabilities, and even if they were, a set of sampled probabilities do not form a probability distribution. Consider if they did - then the first value (itself a sample of size 1) would have to be "1" to make that initial "distribution" integrate to 1, so all subsequent values would have to be "0". You seem to have a very mistaken notion of what's going on, but there are several errors at once, so it's difficult to start untangling your notions. Apr 4, 2015 at 0:23

Quoting from our Bayesian Essentials with R book,

if $\mathscr{D}_n$ denotes a normal $\mathscr{N}\left(\mu,\sigma^{2}\right)$ sample of size $n$, if $\mu$ has a prior equal to a $\mathscr{N}\left(0,\sigma^{2}\right)$ distribution, and $\sigma^{-2}$ an exponential $\mathscr{E}(1)$ distribution, the posterior is given by \begin{align*} \pi((\mu,\sigma^2)|\mathscr{D}_n) &\propto \pi(\sigma^2)\times\pi(\mu|\sigma^2)\times f(\mathscr{D}_n|\mu,\sigma^2)\\ & \propto (\sigma^{-2})^{1/2+2}\, \exp\left\{-(\mu^2 + 2)/2\sigma^2\right\}\\ & \times (\sigma^{-2})^{n/2}\,\exp \left\{-\left(n(\mu-\overline{x})^2 + s^2 \right)/2\sigma^2\right\} \\ &\propto (\sigma^2)^{-(n+5)/2}\exp\left\{-\left[(n+1) (\mu-n\bar x/(n+1))^2+(2+s^2)\right]/2\sigma^2\right\}\\ &\propto (\sigma^2)^{-1/2}\exp\left\{-(n+1)[\mu-n\bar x/(n+1)]^2/2\sigma^2\right\}\,.\\ &\times (\sigma^2)^{-(n+2)/2-1}\exp\left\{-(2+s^2)/2\sigma^2\right\}\,. \end{align*} Therefore, the posterior on $\theta$ can be decomposed as the product of an inverse gamma distribution on $\sigma^2$, $$\mathscr{IG}((n+2)/2,[2+s^2]/2)$$ which is the distribution of the inverse of a gamma $$\mathscr{G}((n+2)/2,[2+s^2]/2)$$ random variable and, conditionally on $\sigma^2$, a normal distribution on $\mu$, $$\mathscr{N} (n\bar x/(n+1),\sigma^2/(n+1)).$$ The marginal posterior in $\mu$ is then a Student's $t$ distribution $$\mu|\mathscr{D}_n \sim \mathscr{T}\left(n+2,n\bar x\big/(n+1),(2+s^2)\big/(n+1)(n+2)\right)\,,$$ with $n+2$ degrees of freedom, a location parameter proportional to $\bar x$ and a scale parameter almost proportional to $s$.

From this distribution, you get the expectation $n\bar x/(n+1)$ that acts as your point estimator of $\mu$. And a credible interval on $\mu$ $$\left(n\bar x/(n+1)-((2+s^2)/(n+1)(n+2))^{1/2}q_{n+2}(\alpha),n\bar x/(n+1)+((2+s^2)/(n+1)(n+2))^{1/2}q_{n+2}(\alpha)\right)$$where $q_{n+2}(\alpha)$ is the $t_{n+1}$ quantile.

• Hey thanks, probably I wasn't clear enough. I understand the derivation of it. It is the computational part that confuses me. I mean if I have say 10 data samples. Now how do I estimate population mean using matlab?
– Sudh
Apr 3, 2015 at 15:44
• You have all the statistical elements above, thus only to replace the three parameters of the Student't density with the values obtained from your sample. If this is a matlab question, it should be asked on Stack Overflow, not here. Apr 3, 2015 at 19:37
• Not specifically matlab but any computational method. How would it be done in practice? For example, matlab can generate random numbers from a t distribution based on degree of freedom provided but how would that correspond to my data? I mean [5 6 8] and [56 54 57] both have 2 degree of freedom but their means are in very different range.
– Sudh
Apr 4, 2015 at 0:25
• Why would you need a random generator? If I have $n=10$ observations with $\bar x=2.3$ and $s^2=312$, the posterior on $\mu$ is a $t(12,2.09,2.38)$ distribution. End of the story. Apr 4, 2015 at 8:41