Mixture coefficients and Parameter estimation for moving Normal distributions

Question

I am searching for a way to do the following steps below. I am stumped and I have stumped a few profs too. Any help would be appreciated.

1. I would like to import data from, say, excel or a .txt file
2. specify a certain number of data points (say 1 to N) to start out with and calculate: a) the coefficients and parameter estimates for a mixture or normal distributions. b) Each mixture component would have the same mean for those data points selected.
3. This process would be carried out now for data points 1 to N+1 until the end of the data set.
4. The values of the mixture coefficients as well as all mean and variance parameters are listed as output.

Essentially I am trying to have a mixture of specific normal distributions estimated a subset of my data one time step at a time.

Are there any gurus out there who know how to do this?

I appreciate any assistance!

Answer 1

Here is how I would do it.

First, import the data using Import. This is straightforward from the documentation, using either Excel or txt files, so for the purposes of this exercise, I'm going to generate the data in Mathematica instead.

myD = MixtureDistribution[{2, 1}, {NormalDistribution[2, 2], 
NormalDistribution[2, 1/2]}];

For reference, here is the PDF.

Plot[PDF[myD, x], {x, -1, 5}]

Here is some fake data generated from that distribution, and an example of a subset being estimated. I understood from your question that you were assuming that the means of the component distributions were equal. It's not hard to change that, of course (just make the second mu into something else like othermu).

pretenddata = RandomVariate[myD, 500];
EstimatedDistribution[Take[pretenddata, 100], 
 MixtureDistribution[{w1, w2}, {NormalDistribution[mu, sig1], 
   NormalDistribution[mu, sig2]}]]

(* MixtureDistribution[{0.703728, 0.296272}, 
 {NormalDistribution[1.99722, 1.96856], NormalDistribution[1.99722, 0.460136]}] *)

A recursive estimation of the kind you want is most easily done as a table. It was quite slow, taking a few minutes on my machine with 100 total sample starting from a subsample of 60.

Table[EstimatedDistribution[Take[pretenddata, i], 
  MixtureDistribution[{w1, w2}, {NormalDistribution[mu, sig1], 
    NormalDistribution[mu, sig2]}]], {i, 400, 500}]

Personally I think that might be overkill, unless you have very small sample sizes - which makes me wonder whether you can justify such a sophisticated distribution. Perhaps you could get just as good results by partitioning the data into subsets.

Now to estimate across separate subsets, you could use Map (/@) to map the estimation function onto subsets created using Partition, like this (where pretenddata has 10,000 points, not 500 as above):

results = EstimatedDistribution[#, 
   MixtureDistribution[{w1, w2}, {NormalDistribution[mu, sig1], 
     NormalDistribution[mu, sig2]}]] & /@ Partition[pretenddata, 1000]

This takes rather less time than the recursive approach, which would be infeasible for such a large data set. Which ever approach you take, though, once you are done, though, you can use standard Part (shorthand [[]]) notation to extract different parts of the results. For example here is a bar chart of the estimated weight on the first component. In the estimated distributions, the component weights are normalised to sum to 1.

BarChart[results[[All, 1, 1]]]

The estimated parameters of the component distributions can be accessed as, for example, results[[All, 2, 1, 1]] - in any case you need a subpart of results[[All,2]].

I do not have the statistical smarts to work out how to do this if you do not know how many component distributions there are in the mixture. Presumably there is some likelihood ratio test you can use to check how many components you need.