We now look at the concept of a prediction interval which is somewhat related to confidence intervals, however, has a different meaning. A prediction interval tells us a predicted range that a single next observation of data is expected to fall within. This differs from a confidence interval

which indicates how confident we are of a particular parameter that we are trying to estimate. For a given distributional model, the bounds of a prediction interval are always wider than those of a confidence interval, as the prediction interval must account for the uncertainty in knowing the population mean, as well as the spread of the data due to variance.

1

The example that we use is for a sequence of data points ????1, ????2, ????3, ..., which come from a normal distribution and are assumed i.i.d. Further assume that we observed ????1,...,???????? but have not yet observed ????????+1. Note that we use ‘little’ ???? for values observed and ‘upper case’ ???? for (yet) unobserved random variables.

In this case, a 100(1 − ????)% prediction interval for the single future observation, ????????+1, is given by,

????̅ − ????1−????,????−1????√1 + 1 ≤ ????????+1 ≤ ????̅ + ????1−????,????−1????√1 + 1 2????2????

where, ????̅ and ???? are respectively the sample mean and sample standard deviation computed from ????1, ... , ????????. Note that as the number of observations, ????, increases, the bounds of the prediction interval decreases towards,

????̅−????1−????????≤????????+1 ≤????̅+????1−???????? 22

To facilitate the task for you and to provide you an idea how you will perform the implementation in R, Python and Octave, the Julia is provided. In the Julia code we illustrate the use of prediction intervals based on a series of observations made from a normal distribution. We start with ???? = 2 observations and calculate the corresponding prediction interval for the next observation. The sample size ???? is then progressively increased, and the prediction interval for each next observation calculated for each subsequent case. The output of the Julia is the Figure 1 which illustrates that as the number of observations increases, the prediction interval width decreases.

1. 1 using Random, Statistics, Distributions, Plots; pyplot()

2. 2 Random.seed!(0)

3

4. 4 mu, sig = 50, 5

5. 5 dist = Normal(mu, sig)

6. 6 alpha = 0.01

7. 7 nMax = 40

8

9 observations = rand(dist,1) 10 piLarray, piUarray = [], [] 11
12 for _ in 2:nMax

13 xNew = rand(dist)
14 push!(observations,xNew) 15

2

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

xbar, sd = mean(observations), std(observations) n = length(observations)
tVal = quantile(TDist(n-1),1-alpha/2)

delta = tVal * sd * sqrt(1+1/n) piL, piU = xbar - delta, xbar + delta

push!(piLarray,piL); push!(piUarray,piU)

end

scatter(1:nMax, observations, c=:blue, msw=0, label="Observations")

plot!(2:nMax, piUarray, c=:red, shape=:xcross, msw=0, label="Prediction Interval")

plot!(2:nMax, piLarray, c=:red, shape=:xcross, msw=0, label="", ylims=(0,100), xlabel="Number of observations", ylabel="Value")

In line 4-7 we setup the distributional parameters, choose, and also set the limiting number of observations we will make, nMax. In line 9 we sample the first sample observation and store it in the array observations. In line 10, we create the arrays piLarray and piUarray, which will be used to store the lower and upper bounds of the prediction interval. Lines 12-23 contain the main logic of this example where in lines 13-14 a new date point is obtained and stored, in lines 16-18 updated prediction interval is calculated and in line 22, the prediction interval is stored for plotting afterwards. Lines 25-30 create Figure 1. Observe the use of shpae=:xcross for setting tick marks.

Figure 1: As the number of observations increases, the width of the prediction interval decreases to a constant.

3

Your task in this is to write programs in Julia, Python, R, and Octave to plot Figure 1. The Julia code is already provided. You are recommended to rely on the concept of this implementation on Julia to reuse the concepts in Python, Octave, and R.

Which programming language that provided relevant solution for this problem? Justify your answer.