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.
 A: 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.
