Search code examples
matlaboctavegaussianmixture-model

What would be valid starting conditions to pass to fitgmdist?


I am trying to use the fitgmdist function from the statistics package in Octave. It works when I let it use the default k++ method for finding initial parameters. However, the results are not consistent and sometimes just plane wrong. That is why I wanted to be able to pass initial values for the means based on looking at the histograms. But for some reason the function won't accept them. The error I keep ending up with is the following:

error: fitgmdist: invalid start parameter
error: called from
    fitgmdist at line 202 column 9
    curve_fitting at line 78 column 17

curve_fitting is just the name of my script. I am trying to use the following code for the function:

nbOrientations = 2;
initial_orientations = [38.0; 18.0]; % #values here should match nbOrientations
initial_weights = ones(1,nbOrientations)/nbOrientations;
initial_Sigma = ones(1,1,nbOrientations);
start = struct('mu',initial_orientations,'Sigma',initial_Sigma,'ComponentProportion',initial_weights)
GMModel_Theta = fitgmdist(Angle_Theta, nbOrientations,'Start', start,'RegularizationValue',0.0001)

My data is just a 700ish by 1 array.

I checked my struct and it seems to me that it satisfies the requirements I could find in the matlab/octave documentation. I am all out of ideas on how to fix this. Hopefully someone can point me in the right direction.

EDIT: I managed to test my script in matlab on someone else's computer and there is just worked. It seems to me that this is an issue with Octave.


Solution

  • Converting the discussion in the comments to an answer for the sake of future readers. There are two issues here:


    The first is that you are encountering a known bug in octave's statistics package v1.4.2, which was reported here, and has been fixed for the upcoming version (which, at the point of writing this answer, has yet to be released).

    If you'd like to apply the fix yourself rather than wait for the next release, effectively correct the typo on line 194 (from 'ComponentProprition' to 'ComponentProportion' ), and comment out the unnecessary check in lines 204-206.


    The second is that you have also come across an unreported bug. I'm converting your code to a full testcase below to demonstrate the issue (I changed the values slightly to match my input):

    pkg load statistics
    
    Angle_Theta          = [ 30 + 10 * randn(1, 10),  60 + 10 * randn(1, 10) ].';
    nbOrientations       = 2;
    initial_orientations = [38.0; 18.0];   % values here should match nbOrientations
    initial_weights      = ones( 1, nbOrientations ) / nbOrientations;
    initial_Sigma        = 10 * ones( 1, 1, nbOrientations );
    
    start = struct( 'mu'                 , initial_orientations,
                    'Sigma'              , initial_Sigma       ,
                    'ComponentProportion', initial_weights        )
    
    GMModel_Theta = fitgmdist( Angle_Theta          , 
                               nbOrientations       ,
                               'Start'              , start ,
                               'RegularizationValue', 0.0001   )
    

    Line 197, tries to ensure that there are no mismatched dimensions. Unfortunately in doing so, it seems to disregard that Sigma may not be shared, so the check fails when sigma contains more than 2 dimensions (i.e. when the 3rd dimension represents the number of components).

    I modified the code, changing size(Sigma) into size(Sigma,1), i.e. effectively making the check against only the rows of Sigma, assuming (naively) that the remaining dimensions of Sigma are fine. This enables the check to pass (while still being a useful check), and the code now runs as expected, giving the following output:

    Gaussian mixture distribution with 2 components in 1 dimension(s)
    Clust 1: weight 0.450954
            Mean: 60.839
            Variance:45.190
    Clust 2: weight 0.549046
            Mean: 32.3048
            Variance:98.217
    AIC=174.207 BIC=179.186 NLogL=82.1037 Iter=84 Cged=1 Reg=0.0001
    

    Since you've helped uncover another bug, it would be useful to report it to the octave bug tracker.

    I'm happy to do it on your behalf; I'm equally happy if you'd like to take this opportunity to engage with the octave / open-source community and contribute the bug report yourself1. Let me know :)


    1. If you do, would you mind commenting a link to the bug report here for reference :)