% MS-A0503 L4A demo: Trying to maximize the likelihood
%
% x = observed data
% family: one of these
%   1 = normal N(mu, 2)     [mu free, sigma fixed]
%   2 = normal N(5, sigma)  [mu fixed, sigma free]
%   3 = uniform Unif(0, b)  [left end fixed, right end free]
% loga: if true, we plot the logarithm of likelihood
%
% The data are plotted as points.
% Then the user can try different parameter values
% by clicking with mouse. End by right click.
%
function paraest(x, family, loga)
    clf
    if nargin<1 || isempty(x)
        % fixed example data
        x = [2 3 6 9];
    end
    if nargin<2
        family = 3;
    end
    if nargin<3
       loga = false;
    end
    
    fprintf('data has average    %.4f\n', mean(x))
    fprintf('data has std.dev.   %.4f\n', std(x,1))
    fprintf('data has maximum    %.4f\n', max(x))
    drawagain(x,family,nan,loga,true,false)

    % this vector collects all parameter values tried
    PARA = [];
    
    while true
        % take next input from user
        [para,~,button] = ginput(1);
        if button ~= 1, break; end
        
        if family==2, para=abs(para); end
        
        drawagain(x,family,para,loga,true,true);
        PARA = [PARA para];
    end
    
    % finally try many more parameter values
    para=linspace(min(PARA), max(PARA), 100);
    drawagain(x,family,para,loga,false,true,true);
end



function drawagain(x,family,para,loga,drawdata,drawlik,connect)
    if nargin<4, loga=false; end
    if nargin<5, drawdata=false; end
    if nargin<6, drawlik=false; end
    if nargin<7, connect=false; end
    
    if nargin>=3
        if family==1
            a = para;
            b = repmat(2,size(para));
        elseif family==2
            a = repmat(5,size(para));
            b = para;
        elseif family==3
            a = repmat(0,size(para));
            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]);

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

            % plot the density function with current parameters
            xgrid = linspace(minx, maxx, 1000);
            if family==3
                fgrid = unifpdf(xgrid,a,b);
            else
                fgrid = normpdf(xgrid,a,b);
            end
            plot(xgrid,fgrid,'-r','linewidth',2);
            
            for xi=x
                if family==3
                    fxi = unifpdf(xi,a,b);
                else
                    fxi = normpdf(xi,a,b);
                end
                plot([xi xi],[0 fxi],'-r');
                plot(xi,fxi,'.r','markersize',20);
            end

            plot([minx maxx],[0 0],'k-');
        end
    end
    
    if drawlik
        % plot the value of the likelihood function
        L = nan(size(a));
        for i=1:length(a)
            if family==3
                L(i) = prod(unifpdf(x, a(i), b(i)));
            else
                L(i) = prod(normpdf(x, a(i), b(i)));
            end
        end
        if loga, L=log(L); end
        title('data and density')
        xlabel('data')
        ylabel('density')
        
        subplot(1,2,2)
        plot(para, L, 'r.', 'markersize', 20);
        hold on
        if connect
            plot(para, L, 'r-', 'linewidth', 2);
        end
        grid on
        if loga
            title('log likelihood');
            xlabel('parameter value');
            ylabel('l');
        else
            title('likelihood');
            xlabel('parameter value');
            ylabel('L');
        end

        if connect
            % We did a grid search so report the results
            [~, imax] = max(L);
            fprintf('best from grid search: L(%.4f) = %.9f\n', ...
                para(imax), L(imax))
        else
            % We only plotted one attempt
            fprintf('L(%.4f) = %.9f\n', para, L);        
        end
    end
end