Mixture model

Normal distributions mixture model/ Waiting times of eruptions of the Old Faithful

This loads the package.

<<StatisticalInference`

The following table gives the waiting time (in minutes) of 299 consecutive eruptions of the Old Faithful  geyser in Yellowstone National Park (Pawitan 2001, p. 350).

dataFile = ToFileName[
{$InstallationDirectory,
"\\AddOns\\Applications\\StatisticalInference\\Datasets"},
"Eruptions.dat"];
dataOnEruptionWaitingTimes = Import[dataFile,"List"]

This is a list plot of waiting times.

ListPlot[dataOnEruptionWaitingTimes,Axes -> False,Frame -> True,PlotRange -> All];

In the list plot there is an indication of a gap in the vicinity of the waiting time 65, which is confirmed by the empirical cumulative distribution function.

Plot[Evaluate[DistributionFunction[EmpiricalModel[dataOnEruptionWaitingTimes],y,{}]],{y, 40, 110}];

Response consists of the waiting times of eruptions.

resp=dataOnEruptionWaitingTimes;

Assume that the waiting times can be modelled as a sample from a mixture of two normal distributions.

ℳ = 
        SamplingModel[ 
          MixtureModel[
            NormalModel[μ, σ], 2], Length[resp]];

This calculates the fitted model.

ℱ = 
      MLEFit[ℳ, resp, 
        StartingPoint -> {{0.3, 0.7}, {{55, 4}, {80, 8}}}]

The maximum likelihod estimate is

MatrixForm /@ MLParameterEstimate[ℱ]  // N

This shows the fitted and empirical cumulative distribution functions in same graph.

Plot[Evaluate[{DistributionFunction[EmpiricalModel[dataOnEruptionWaitingTimes],y,{}],DistributionFunction[
                  MixtureModel[NormalModel[μ, σ], 2], y, 
                  MLParameterEstimate[ℱ]]}],{y, 40, 110}];

This calculates (approximate) 95%-level profile likelihood based confidence interval for the probability of the first component of the mixture.

ℐ = 
      ProfileInterval[ℱ, 
        p, {{p, 1 - 
              p}, {{μ1, σ1}, {μ2, σ2}}}]

ConfidenceLimits[ℐ]

This calculates (approximate) 95%-level profile likelihood based confidence interval for the probability of the ratio of variances.

ℐ = 
      ProfileInterval[ℱ, σ((1^2)/σ)(2^2), {{p, 
            1 - p}, {{μ1, σ1}, {μ2, σ2}}}]

ConfidenceLimits[ℐ]

This calculates the likelihood ratio test for the hypothesis of equal variances.

𝒯 = 
      LRTest[ℱ, {p, μ1, μ2, σ}, {{p, 
            1 - p}, {{μ1, σ}, {μ2, σ}}}, 
        StartingPoint -> {0.3, 64, 80, 6}]

The observed significance level is

ObservedSignificanceLevel[𝒯]