Sannsynlighetsmaksimering

Likelihoodfunksjon for normalfordelte data med kjent varians. Normaltettheten kan flyttes fram og tilbake (dvs. at µ varieres), og høydene fra datapunktene på x-aksen og opp til tettheten vises. Likelihoodfunksjonen er produktet av disse høydene. Fungerer for MATLAB R2016a.

NB! Filen mle.m (nederst) må ligge i katalog spesifisert i den øverste koden (linje som begynner med cd).

% n uavhengige høyder, som er N(my,6^2) skal brukes til å estimere my = 180 ved
% sannsynlighetsmaksimering. Likelihoodfunksjonen er produktet av høydene over hvert
% datapunkt opp til tettheten til N(my,6^2).
% 
% Dra "slideren" med mus (eller pilknapper hvis slideren har fokus) for å variere my og
% vise høydene og likelihoodfunksjonen.
%
% Kopier alt dette inn i kommanduvinduet til MATLAB, og pass på at filen mle.m ligger i
% katalogen spesifisert i neste linje.
%
cd('n:/bakke/pc/Undervisning/MATLAB'); % katalog for mle.m
global n;
n=30; % kan varieres
global x;
x=normrnd(180,6,1,n);
global sigma;
sigma=6;
global lnlgrense;
lnlgrense=-100;
clf;
slider=uicontrol('Style','slider','Min',160,'Max',200,'Units','normalized',...
 'Position',[.04 .00 .82 .05],'Callback',@mle);
% For å få passende skala for lnL:
set(slider,'Value',200);
ekstrem2=mle(slider);
set(slider,'Value',160);
ekstrem1=mle(slider);
ekstrem=-min(ekstrem1,ekstrem2);
k=floor(log10(ekstrem));
if floor(ekstrem*10^(-k)+1)==1
  k=k-1;
end
lnlgrense=-10^k*floor(ekstrem*10^(-k)+1);
% Så startbildet:
mle(slider);

Filen mle.m:

% Dette er filen mle.m
function [lnl] = mle(source,callbackdata)
  global n;
  global x;
  global sigma;
  global lnlgrense;
  my=get(source,'Value');
  subplot(1,10,1:9);
  plot(160:.1:200,normpdf(160:.1:200,my,sigma),'color',[1 0 0]);
  xlim([160 200]);
  set(gca,'Position',[.1 .13 .6956 .8150],'Units','normalized');
  hold on;
  scatter(x,repmat(0,1,n),'filled','MarkerFaceColor',[.5 .5 1],'MarkerEdgeColor',[0 0 1]);
  tekst2=('m');
  title(['\mu = ' num2str(my)]);
  hold off;
  lnl=0;
  for i=1:n
    line([x(i) x(i)],[0 normpdf(x(i),my,sigma)])
  end;
  lnl=-.5*sum(((x-my)/sigma).^2)-n*(.5*log(2*pi*sigma^2));
  subplot(1,10,10);
  bar(1,lnl,1);
  ylim([lnlgrense 0])
  set(gca,'XTickLabel','','XTick',0,'Position',[0.9 0.13 0.03 0.8150],'Units','normalized');
  handle=title(['ln\it L\rm\bf = ' num2str(lnl)]);
  set(handle,'Position',[0 .5629 0]);
  hold off;
end
2017-03-16, Øyvind Bakke