function testParameterEstimation3
%% prepare sample sizes and true parameter values
k=4;
nX=5;
nY=7;
meX=2*rand(1,k)-1;
meY=2*rand(1,k)-1;
cov0=cov(2*rand(k+1,k));
%% calculate true discriminant function parameters
V=ones(1,k+1);
V(1,2:k+1)=(meY-meX)/cov0;
V(1,1)=V(1,2:k+1)*(meY+meX)'/2;
%% Monte-Carlo
L=1000;
res=ones(L,k+1);
for i=1:L
res(i,:)=func(k,nX,nY,meX,meY,cov0,V);
end;
%% plot results
for j=1:k+1
figure(j);
hist(res(:,j),50);
hold on;
[nn,xx]=hist(res(:,j),50);
textt=['parameter ',num2str(j),newline];
textt=[textt,'true value ',num2str(V(j)),newline];
textt=[textt,'mean value ',num2str(mean(res(:,j))),newline];
textt=[textt,'skewness ',num2str(skewness(res(:,j)))];
text(xx(5),max(nn)/2,textt);
hold off;
end;
function r=func(k,nX,nY,meX,meY,cov0,param0)
% generate X and Y, calculate main parameter estimation
r=ones(1,k+1);
X=mvnrnd(meX(ones(nX,1),:),cov0);
Y=mvnrnd(meY(ones(nY,1),:),cov0);
mX=mean(X);
mY=mean(Y);
covX=cov(X);
covY=cov(Y);
S=((nX-1)*covX+(nY-1)*covY)/(nX+nY-2);
% calculate discriminant function parameters estimations
Vc=ones(1,k+1);
Vc(1,2:k+1)=(mY-mX)/S;
Vc(1,1)=Vc(1,2:k+1)*(mY+mX)'/2;
r=Vc;