# With a small sample from a normal distribution, do you simulate using a t distribution?

I want to simulate temperature data for some "what-if" calculations. The problem is, I only have a time series of 10 actual temperature data values. I want to use temperature as an input to the simulation, so I need a way to generate a large number of temperature values that are consistent with the original 10 values. It's probably ok to assume that they came from a normal distribution, but I don't know the mean or the variance.

I have no way to prove it, but I doubt the 10 values do a good job of representing the full temperature range. If I use the sample function for the simulation, as shown below, I only get the original values back. That just doesn't look right. If I use the rnorm function, I know that I don't know the variance, so I don't think that is right either. So, I'm left with the rt function (t-distribution).

Below is a mock up of the problem.

ori <- rnorm(n=10, mean=65, sd=5) #original 10 data points

num.sam <- sample(x=ori, size=100, replace = TRUE) #simulation using sample
num.tdis <- mean(ori) + (rt(n=100, df=10) * sd(ori)) #simulation using a t distribution

hist(ori, breaks=40:90)
hist(num.sam, breaks=40:90)
hist(num.tdis, breaks=40:90)


My questions are,

1. When I only have data (mean and variance unknown), and it is reasonable to assume that the data came from a normal distribution, is it ok to generate data for a simulation using a t-distribution?

2. For this type of situation, the only time I would use rnorm for the simulation is if I knew the variance (not a variance estimated from the data), right?

3. If a t-distribution simulation is ok for these conditions, are there any conditions where it is better to just sample the data (for example 100 original data points, 200, etc)?

Edit:

1. Since I used the original data to estimate the mean and variance, should the degrees of freedom in the third line of the code (for rt(...)) be reduced from 10 to 9? Or 8?
• What's the purpose of the simulation? – mark999 Nov 18 '12 at 20:47
• I'll use the simulated temperature as an input to feed other calculations, so I can estimate a spread for the final forecast. I have a related issue about a moving temperature mean and possibly a moving variance, that I'll ask as another question, but this is a simplified version of the main issue. – xls5929 Nov 18 '12 at 22:24
• I still don't understand what you want to do, but maybe you could use simulations to compare what happens if you simulate using a t distribution compared with using a normal distribution. – mark999 Nov 19 '12 at 9:05
• I'm simply trying to generate a large number of temperature values that are consistent with the original 10 value sample. The reason that I assembled the above mock-up code was to compare the generated rt() values (not shown, but I also looked at rnorm() values) with the ori values. In the above mock-up, I actually know the mean and variance of ori. The problem was/is, nothing stood out as the "right" way to simulate the data. That's why I'm asking if someone can explain what is the "right" way and why. I may be wrong in assuming this is a common question. – xls5929 Nov 19 '12 at 15:59

Using your assumptions that the points came from a normal distribution with unknown mean and variance, The T distribution is the correct distribution to sample from no matter how many data points you have because it is the posterior predictive distribution of your model. You might want to check your formula though as it looks a bit simpler than what I've seen before.

• I tested/re-tested/tested-again the formula (line 3, where rt() is used) and I'm not sure what the problem is. Can you be more specific? Are you talking about a different line in the code? I've been digging around on the net, and apparently simulating with a t-distribution is not as common as I originally expected. – xls5929 Nov 21 '12 at 20:44
• See here for example and note the posterior predictive distribution. You are not "simulating", but rather precisely following your Gaussian modelling assumption. There's a nice derivation here. – Neil G Nov 22 '12 at 9:44

You could generate a vector of means from a normal distribution (or t if you prefer) representing your uncertainty in the mean, then generate a vector of variances from a $\chi^2$ distribution representing your uncertainty in the variance, then generate the actual observations from a normal with your vector of means and the vector of variances as the parameters. This will take into account the extra levels of uncertainty that you mention.

If you have some feel for where you think the mean and/or variance should be (but don't know exactly) then you may want to try a Bayesian approach where you can use that prior information.

• I tried your method, but the spread was bigger than my method above, and I wasn't real sure that I did it right. I generated random sigmas using the following code sigma <- sqrt((datpoi - 1) * (sd(ori)^2)/rchisq(n=newdatpoi, df=datpoi-1)) where datpoi is the number of original data points (in my example, datpoi=10). And, newdatpoi is the number of random values wanted (in my example, newdatpoi=100). Those sigmas and some random means (centered around mean(ori) ) were fed into rnorm() to generate the new data points. As far as the Bayes approach goes, that's out of my league. – xls5929 Nov 21 '12 at 20:56

I would elaborate Neil G and Greg Snow's answers as follows :

• run a noninformative Bayesian inference for your original $10$ data values
• use the posterior predictive distribution to generate new data

The posterior predictive distribution derived from a noninformative prior exactly aims to provide your desire: a distribution that generates data "consistent with the original data", taking into account the uncertainty about the model parameters.

Now, what is the posterior predictive distribution derived from the noninformative prior ? That depends on the choice of the noninformative prior, but there is a good "default" noninformative prior for the normal sample model. You can also "cheat" a little and use the "Bayesian-frequentist" predictive distribution (also called sometimes "the frequentist predictive distribution"). The principle of the frequentist predictive distribution is the following one. The classical $100(1-\alpha)\%$-prediction interval for a new observation is $\bar{y} \pm \mathrm{t}^*_{n-1}(\alpha/2) \hat\sigma\sqrt{1+\frac{1}{n}}$. Then the Bayesian-frequentist predictive distribution is taken to be the distribution of $\bar{y} + T \hat\sigma\sqrt{1+\frac{1}{n}}$ where $\bar{y}$ and $\hat\sigma$ are considered as fixed and $T$ has the Student $\mathrm{t}_{n-1}$ distribution. Thus, the $100(1-\alpha)\%$-quantile of the frequentist predictive distribution equals the usual $100(1-\alpha)\%$-upper prediction bound.

I do not exactly remember the Bayesian predictive distribution derived from the default noninformative prior but it is very close to the frequentist predictive distribution (there are some slight differences such as $\mathrm{t}^*_{n-\frac{1}{2}}$ instead of $\mathrm{t}^*_{n-1}$). I will update my answer when I will find the formulas.

Here I asked a question related to the performance of these predictive distributions.

I claimed that the frequentist predictive distribution is derived from "little cheating" because it does not really has a theoretical fundation. But I'm sure it is possible to show the performance of the use of this distribution in a frequentist sense.