I am rather new to Matlab and never had a lot to do with statistics, so I apologize already for possibly being ignorant of quite a bit of important knowledge. It also would be nice if you could answer as simple as possible. Thanks ☺

My question is about finding the parameters of a univariate distribution. I get two sets of data, which I want to analyze automatically. The first set is bimodal, with the two distributions generally being separated completely. The second set has at least three modes. My data is about the fiber directions in a carbon composite. So if I have 2 general directions I get the bimodal distribution, but if I have 3, I get two high peaks and usually one or two smaller ones (the peak at 0° is wrong data I'm working on eliminating).

I already played around a bit with mle, and got it to work for the bimodal case, but it takes ages, even when given good starting values. Right now I am splitting the bimodal in two and fit two separate normal distributions, but I would like to find a more elegant way. As for the 3-fiber case, I’m still without any working idea.

I really appreciate any help with this :)

The distributions

enter image description here

enter image description here

  • $\begingroup$ Just to be clear at this. You don't just want the functionality of fitdist. Do you? Also to make sure about this : You are trying to find the location and the amplitude of the distribution's modes, if I understand? Do you use the ksdensity function? $\endgroup$ – usεr11852 says Reinstate Monic Jun 18 '13 at 11:19
  • 1
    $\begingroup$ On second thought: Won't gmdistribution.fit(data,k);, $k$ being the number of peaks you expect to have in this case, give you what you want? $\endgroup$ – usεr11852 says Reinstate Monic Jun 18 '13 at 11:52
  • $\begingroup$ Hello user11852, thanks for your quick answer. I was actually expecting to get an email alert, so I didn't check sooner, thus the late reply. Sorry! I actually found the ksdensity function and it works well smoothing my data, and I get my peaks through it. From there I work with different cases and try fitting different distributions on parts of the data. I would though like to find a better solution. $\endgroup$ – Elise Jun 26 '13 at 10:21
  • $\begingroup$ As for the gmdistribution: It sounds interesting, but I don't think I really understand it. At first, fitting it on the bimodal data, I got a 2 component, 703 dimension obj? And after that it always gave back an error, saying something about Ill-conditioned covariance. Do I maybe need to alter my data first? But it worked once? Thanks again! $\endgroup$ – Elise Jun 26 '13 at 10:24
  • $\begingroup$ Ill-conditioned covariance will probably refers to collinearity effects in your data but I can't be sure. Aren't your data a 1-D vector? I do not think you need to alter you data; "maybe" centre them to zero. Could please give an example where this happens (possibly by editing your post?) Is it just the trimodal case that is the problem and not the bimodal? What part of the gmdistribution.fit(data,k) confuses you? $\endgroup$ – usεr11852 says Reinstate Monic Jun 26 '13 at 10:42

Just to give a worked example of what I am talking about in simple case of a trimodal distribution:

rng( 0 ,'twister')  %Set seed to 0
%Set your starting hyperparameters
Modality = 3;
Sigmas = [ 2 1 3]; %Standard Deviations
Means =  [-5 1 9];
MixCoefs = [.4 .25 .35]; 

Length = 1000;
%Get matrix with individual distributions (non-effective way just for illustr.)
X = randn(Length,Modality).*repmat(Sigmas,Length,1)+repmat(Means,Length,1);
%Mix the distributions to your MixCoefs to get the final mixture.
Y = [ X(1: (Length* MixCoefs(1)),1); 
      X(1: (Length* MixCoefs(2)),2); 
      X(1: (Length* MixCoefs(3)),3) ]; 

%Uncomment to visually inspect your empirical pdf and check multimodality

%Fit a trimodal to your data
fitted_obj3 = gmdistribution.fit(Y, 3);
%Fit a bimodal to your data
fitted_obj2 = gmdistribution.fit(Y, 2);

%Check which if the bimodal fit is better based on AIC
fitted_obj2.AIC < fitted_obj3.AIC
%Check which if the bimodal fit is better based on BIC
fitted_obj2.BIC < fitted_obj3.BIC
%Unsurprisingly the trimodal is better is both cases.

best_fitted_obj = fitted_obj3; %Watch it the ordering probably it is not 
% the one you started with

Est_Sigmas = sqrt( best_fitted_obj.Sigma);  %[ 1.03  3.00  2.00]
Est_Means =  best_fitted_obj.mu;            %[ 1.09  9.15 -5.11]

%You can get the  mixing proportions directly by using:

%Or get an estimate for each reading about the possibility of being part in
%a specific distribution by using the .posterior() functionality :

%Posterior "ownership" probabilities of a "-11" reading 
%ans =
%    0.0000    0.0000    1.0000 %So pretty certainly it is on the third
%    component (unsurprisingly) so from the estimated N( -5.11 , 2.00^2)

%Posterior "ownership" probabilities of a "3" reading 
%ans =
%    0.7557    0.2434    0.0009 %Probably on the first component ie. from
%    the N( 1.09, 1.03^2) but actually the second distro is not too
%    unlikely either

Clearly more readings will give you a better fit and less readings a worse one (on average at least). I tried to extensively comment the code so it is easier to follow (as a general coding advice don't comment absolutely everything cause comments require maintenance also...). Computationally speaking gmdistribution.fit is using an E-M algorithm to find the solution but I don't think that is an issue for you. Hope this code helps a bit. If you have questions fire away. :)

  • $\begingroup$ You're amazing. Thank you so much for writing this up for me! It actually helped a quite a bit in my understanding, I think, also some things are still vague. So the questions: Why did you use the square root to get Est_Sigma? And what are these posterior probabilities? I read the matlab article, but I don't think I understood what it was about. Is it the probability of where data can lay regarding each component? $\endgroup$ – Elise Jun 26 '13 at 13:38
  • $\begingroup$ 1. the .Sigma stores the variances, you want the standard deviations. 2. Effectively you are correct. It returns a vector of length $k$ that quantifies the probability that a point $x$ is part of component $i$, $i = \{1,2..,k\}$ given you have $k$ components. That's why here we get 3 numbers adding up to 1; any given reading must be part of a component. $\endgroup$ – usεr11852 says Reinstate Monic Jun 26 '13 at 14:05

Not the answer you're looking for? Browse other questions tagged or ask your own question.