gmdistribution.fit and gmdistribution help needed

I tried to model a multivariate gaussian density with just a data set to estimate the mean, covariance and mixing parameter using gmdistribution.fit. But i dont know whether its correct. Here is my code:

function Ecc = Econtrol(O,K,m,T,n,q1,p2)
x = reshape(O(1:2*n),2,n);
U1 = reshape(O(2*n+1:q1),1,p2);
x=x';
obju = gmdistribution.fit(U1',K,'SharedCov',true,'CovType','diagonal');
objx = gmdistribution.fit(x,K,'SharedCov',true,'CovType','diagonal');

px=0;
for k=1:m
px = log(pdf(objx,x(k,:))+pdf(objx,x(k,:)))+px;
end
pu=0;
for k=1:T-m
pu = log(pdf(obju,U1(:,k))+pdf(obju,U1(:,k)))+pu;
end
Ecc = -px -pu;
end

below is the equation i wanna model. is it correct?

%

[Daniel Shub]

Closed as doit4me, please show your what you have tried and where you are stuck.

Answer 1

[Tom Lane]

I don't know if your code is correct, but:

1. The text seems not to declare that the covariance is diagonal, and does use notation to suggested that it may not be shared across mixture components.

2. Unless I'm just not seeing something, you seem to compute pdf(objx,x(k,:)) twice inside the log, and the same for U, but I don't know why you would want that.

3. Your text describes U as multivariate, but you seem to create U1 with just one row and index into it repeatedly, getting scalar values each time.

4. It may be possible to avoid the loops and compute the pdf on an entire array in one call, then sum the log of the result.

Answer 2

[Wei Cai Law]

Hi Tom, thanks for your reply. I am still new in matlab hence i am making quite a few silly mistake.

I tried computing this function without declaring the covariance diagonal and share covariance but an error occurred. This is the error: Ill-conditioned covariance created at iteration 2. Is there a way to solve this error?

In this multivariate model, i have two component hence i computed pdf(objx,x(k,:)) and for U twice. I tried editing it. Below is the edited version. Does it makes more sense?

obj_u = gmdistribution.fit(U1',K);
obj_x = gmdistribution.fit(x,K);
mu_u = obj_u.mu;
mix_u = obj_u.PComponents;
cova_u = obj_u.Sigma;
mu_x = obj_x.mu;
mix_x = obj_x.PComponents;
cova_x = obj_x.Sigma;

px=0;
for k=1:m
%px = log(pdf(objx,x(k,:))+pdf(objx,x(k,:)))+px;
px = log(mix_x(1)*(mvnpdf(x(k,:),mu_x(1),cova_x(1)))+mix_x(2)*(mvnpdf(x(k,:),mu_x(2),cova_x(2))))+px;
end
pu=0;
for k=1:T-m
%pu = log(pdf(obju,U1(:,k))+pdf(obju,U1(:,k)))+pu;
pu = log(mix_u(1)*(mvnpdf(U1(k,:),mu_u(1),cova_u(1)))+mix_u(2)*(mvnpdf(U1(k,:),mu_u(2),cova_u(2))))+pu;

[Tom Lane]

An ill-conditioned matrix can result from various things. One might be a cluster with too few components to estimate a full-rank covariance. Sharing or forcing a diagonal would be possible ways to deal with that.

pdf(obj,x), where obj is a gmdistribution object, computes the full mixture distribution for you. So I think you want to revert back to the old code, but compute log(pdf()) rather than log(pdf()+pdf()). Furthermore, you don't need to loop over a single row at a time. Check this out, showing how you can compute the pdf for the entire array and get the same answer as if you did each row separately:

>> x = [0 0;0 1;1 1;10 10;10 11;11 11];
>> g = gmdistribution.fit(x,2);
>> pdf(g,[.4 .5]) + pdf(g,[.7 .7])
ans =
    0.6487
>> pdf(g,[.4 .5;.7 .7])
ans =
    0.3631
    0.2856
>> sum(ans)
ans =
    0.6487

[Wei Cai Law]

Hi, thank you for your answer. i find it very useful. can i ask you what are the values in the pdf function, pdf(g,[.4 .5]). so can i say that with pdf(g,[.4 .5;.7 .7]), i can compute the two components? and i would need to increase the number of component to prevent the same error?

[Tom Lane]

I don't understand. The variable g represents both components. pdf(g,x) compute the sum over both. Are you asking how to get at each one individually?

[Wei Cai Law]

i meant how you get the values .4 and .5? and for the ill conditioned matrix, should i increase the number of component to solve the error?