function gamma = func_theoreticalContrast(Q,Mth,O)
% Computes L^2 distance to projection of the real density (i.e. the variance).
% Q:   k*k   transition matrix
% Mth: m*m*m theorical joint distribution
% O:   m*k   coeff of emission laws
%
% Q and O are used to compute the joint distribution Mhat, from which we deduce the variance ||Mhat-Mth||^2.

maux=size(Mth,1);
kaux=size(Q,1);
gamma=0;
pi = func_stationaryLaw(Q);

for b=1:maux
    Mhatb = O * diag(pi) * Q * diag(O(b,:)) * Q * O';
    
    gamma = gamma + sum(sum( (Mhatb - squeeze(Mth(:,b,:))) .* (Mhatb - squeeze(Mth(:,b,:))) ));
end
