test_em_upm.m

clear,clc,close all hidden

disp('A simple example');

d                     = 2;
m                     = 1;
Ntrain                = 1000;
Ntest                 = 1000;
nbite                 = 100;

R                     = 15*rand(m , m , d);
R                     = sort(R);

P                     = rand(m , m , d );
sumP                  = sum(P , 3);
P                     = P./sumP(: , : , ones(d , 1));

[Ztrain , Xtrain]     = sample_upm(Ntrain , R , P);
[Ztest  , Xtest]      = sample_upm(Ntest , R , P);



R0                    = 10*rand(m , m , d);

P0                    = rand(m , m , d);
sumP                  = sum(P0 , 3);
P0                    = P0./sumP(: , : , ones(d , 1));

[logl , Rest , Pest]  = em_upm(Ztrain , R0 , P0 , nbite);
[Rest , ind]          = sort(Rest);
Pest                  = Pest(ind);


Ltrain                = denspoiss(Ztrain , Rest , Pest);
[val , Xtrain_est]    = max(Ltrain , [] , 1);
Err_train             = sum(Xtrain_est ~= Xtrain)/Ntrain;
disp(sprintf('Train error = %f', Err_train));

Ltest                 = denspoiss(Ztest , Rest , Pest);
[val , Xtest_est]     = max(Ltest , [] , 1);
Err_test              = sum(Xtest_est ~= Xtest)/Ntest;
disp(sprintf('Test error = %f', Err_test));

mZ                    = max(Ztrain);
x                     = (0:mZ);
bin                   = histc(Ztrain , x)/Ntrain;
pdf                   = mupoispdf(x , R , P);
pdfest                = mupoispdf(x , Rest , Pest);


figure(1), plot(Ztrain)
figure(2) , bar(x , bin);
hold on
h = plot(x , pdf , 'r' , x , pdfest , 'g' ,  'linewidth' , 2);
hold off
legend(h , 'True' , 'Estimated' , 'location', 'best');
axis([ -2 , mZ + 2 , 0 , 1.2*max(pdf)])