Contents
margpoly.m
From A First Course in Machine Learning, Chapter 2. Simon Rogers, 01/11/11 [simon.rogers@glasgow.ac.uk] Marginal likelihood for model selection
clear all;close all;
Generate the data
N = 100; x = sort(rand(N,1)*10-5); noise_var = 150; t = 5*x.^3 - x.^2 + x; % Try adding and removing terms from this function, or changing term % weights. e.g. % t = 0.0005*x.^3 - x.^2 + x; t = t + randn(size(x)).*sqrt(noise_var); % Plot the data plot(x,t,'k.','markersize',10); xlabel('x'); ylabel('t');
Fit models of various orders
orders = [0:8];
testx = [-5:0.01:5]';
X = [];
testX = [];
for i = 1:length(orders)
si0 = eye(orders(i)+1);
mu0 = repmat(0,orders(i)+1,1);
X = [X x.^orders(i)];
testX = [testX testx.^orders(i)];
siw = inv((1/noise_var)*X'*X + inv(si0));
muw = siw*((1/noise_var)*X'*t + inv(si0)*mu0);
% Plot the data and mean function
figure(1);hold off;
plot(x,t,'k.','markersize',10);
xlabel('x');
ylabel('t');
hold on
plot(testx,testX*muw,'r');
ti = sprintf('Model order %g',orders(i));
title(ti);
% Compute the marginal likelihood
margcov = noise_var*eye(N) + X*si0*X';
margmu = X*mu0;
D = length(margmu);
log_marg(i) = -(D/2)*log(2*pi) - 0.5*log(det(margcov));
log_marg(i) = log_marg(i) - 0.5*(t-margmu)'*inv(margcov)*(t-margmu);
end
Plot the marginal likelihoods
figure(1); hold off bar(orders,exp(log_marg)); xlabel('Model Order'); ylabel('Marginal likelihood');
