Center of gaussian mixed distribution area

Hi! I want to calculate the center of a gaussian mixed distribution. The code goes as follows:

dd = importdata('gmdata.txt'); %%%% import data
xx = dd(:,1);
yy = dd(:,2);

obj = gmdistribution.fit([xx,yy],2); %%%% fit a distribution
figure;
scatter(xx,yy,10,'.') %%%% Scatter plot with points of size 10
hold on
gmPDF = @(x,y) arrayfun(@(x0,y0) pdf(obj,[x0 y0]),x,y);
hh = fcontour(gmPDF,[-8 6]);

I am getting the above distribution as output. I want to recreate the value of GM distribution centers as given in "Tanabe, Hiroko, Keisuke Fujii, and Motoki Kouzaki. "Intermittent muscle activity in the feedback loop of postural control system during natural quiet standing." Scientific reports 7.1 (2017): 10631." The sentences read as:

Kindly help me with the code. I guess the center they are mentioning as the center of the fcontour rings, but not sure.

[Torsten]

Your data ressemble a usual multivariate Gaussian, don't they ?

https://de.mathworks.com/help/stats/mvnpdf.html

Answer 1

[the cyclist]

I downloaded and looked over the paper. I am also not certain what they mean by the "center" of the mixed gaussian model.

But, I think a pretty reasonable guess is they mean the peak of the joint PDF. (I think is also what you mean by the "center of the contour rings").

You can use the fminsearch function to find this peak:

dd = importdata('gmdata.txt'); %%%% import data
xx = dd(:,1);
yy = dd(:,2);

obj = gmdistribution.fit([xx,yy],2); %%%% fit a distribution
figure;
scatter(xx,yy,10,'.') %%%% Scatter plot with points of size 10
hold on
gmPDF = @(x,y) arrayfun(@(x0,y0) pdf(obj,[x0 y0]),x,y);
hh = fcontour(gmPDF,[-8 6]);

% Define the function to be minimized (using the negative because we actually want the maximum)
f = @(x) -obj.pdf(x);


% Initial guess as to where the maximum is.
% (This is a relatively poor guess, but the result looks OK anyway)
xy_init = [0 0]; 

% Use fminsearch to find the maximum
xy_max = fminsearch(f, xy_init)

xy_max = 1×2

-1.9393 -5.2685

Eye-balling the plot, this looks accurate.

[Torsten]

Yes, it's one of the obj.mu (where both of them look pretty much the same).

[Aishwarya Rao]

Thank you @the cyclist. The code worked!