% demonstraatio parametrin estimoinnista
% x = havaittu data
% perhe: valitsee jakaumaperheen näistä
%   1 = normaali N(mu, 2)     [mu vapaa, sigma fiksattu]
%   2 = normaali N(5, sigma)  [mu fiksattu, sigma vapaa]
%   3 = tasajak  Tas(0, b)    [alkupää fiksattu]
% loga: true valitsee logaritmisen uskottavuusfunktion
%
% Kun data on piirretty pisteparvena,
% käyttäjä voi kokeilla eri parametriarvoja hiirellä
% klikkaamalla.
%

%% GNU Octavessa muista ajaa "pkg load statistics" ennen paraest-funktiota.
%% Tämä ei toimi gnuplot-piirturilla.
%% Kokeile "graphics_toolkit('fltk')" tai "graphics_toolkit('qt')" jos ei muuten toimi.

function paraest(x, perhe, loga)
    clf
    if nargin<1 || isempty(x)
        % esimerkkidata
        x = [2 3 6 9];
    end
    if nargin<2
        perhe = 1;
    end
    if nargin<3
       loga = false;
    end
    
    fprintf('datan keskiarvo on    %.4f\n', mean(x))
    fprintf('datan keskihajonta on %.4f\n', std(x,1))
    fprintf('datan maksimi on      %.4f\n', max(x))
    drawagain(x,perhe,nan,loga,true,false)

    % tähän kerätään kaikki inputit
    PARA = [];
    
    while true
        % seuraava kokeiltava mu-arvo käyttäjältä
        [para,~,button] = ginput(1);
        if button ~= 1, break; end
        
        if perhe==2, para=abs(para); end
        
        drawagain(x,perhe,para,loga,true,true);
        PARA = [PARA para];
    end
    
    % lopuksi kokeillaan tiheästi lisää parametriarvoja
    for para=linspace(min(PARA), max(PARA), 100)
        drawagain(x,perhe,para,loga,false,true);
    end
end



function drawagain(x,perhe,para,loga,drawdata,drawlik)
    if nargin<4, loga=false; end
    if nargin<5, drawdata=false; end
    if nargin<6, drawlik=false; end
    
    if nargin>=3
        if perhe==1
            a = para;
            b = 2;
        elseif perhe==2
            a = 5;
            b = para;
        elseif perhe==3
            a = 0;
            b = para;
        end
    end
    
    if drawdata
        subplot(1,2,1)
        cla
        plot(x,0,'b.','markersize',20);
        grid on
        hold on

        minx = -5;
        maxx = 20;
        set(gca,'xlim',[minx maxx]);
        set(gca,'ylim',[0 1]);

        if drawlik
            if perhe~=2
                plot(para,0,'rx','markersize',20);
            end

            % piirrä tiheysfunktio näillä parametreilla
            xgrid = linspace(minx, maxx, 1000);
            if perhe==3
                fgrid = unifpdf(xgrid,a,b);
            else
                fgrid = normpdf(xgrid,a,b);
            end
            plot(xgrid,fgrid,'-r');
            
            for xi=x
                if perhe==3
                    fxi = unifpdf(xi,a,b);
                else
                    fxi = normpdf(xi,a,b);
                end
                plot([xi xi],[0 fxi],':k');
                plot(xi,fxi,'.r','markersize',20);
            end

            % lukusuora
            plot([minx maxx],[0 0],'k-');
        end
    end
    
    if drawlik
        % laske ja piirrä parametrin uskottavuusfunktiota
        if perhe==3
            ff = unifpdf(x, a, b);
        else
            ff = normpdf(x, a, b);
        end
        L = prod(ff);
        if loga, L=log(L); end
        
        subplot(1,2,2)
        plot(para, L, 'r.', 'markersize',20);
        hold on
        grid on
        if loga
            title('log. uskottavuusfunktio');
            xlabel('parametrin arvo');
            ylabel('l()');
        else
            title('uskottavuusfunktio');
            xlabel('parametrin arvo');
            ylabel('L()');
        end
        %fprintf('L(%.4f) = %.9f\n', para, L);        
    end
end