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.

  • $\begingroup$ 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. $\endgroup$ – Xi'an Apr 3 '15 at 7:18
  • $\begingroup$ Why should they add to 1? $\endgroup$ – Glen_b -Reinstate Monica Apr 3 '15 at 8:33
  • $\begingroup$ Shouldn't a probability distribution integrate to 1? $\endgroup$ – Sudh Apr 3 '15 at 15:41
  • $\begingroup$ 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. $\endgroup$ – Glen_b -Reinstate Monica Apr 4 '15 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.

  • $\begingroup$ 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? $\endgroup$ – Sudh Apr 3 '15 at 15:44
  • $\begingroup$ 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. $\endgroup$ – Xi'an Apr 3 '15 at 19:37
  • $\begingroup$ 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. $\endgroup$ – Sudh Apr 4 '15 at 0:25
  • $\begingroup$ 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. $\endgroup$ – Xi'an Apr 4 '15 at 8:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.