1 function sheet08
   2 
   3 % The parameters:
   4 M = 1000; % number of evaluation points
   5 X = linspace(-4, 4, M)'; % the evaluation points
   6 ITERS = 100; % number of resamples 
   7 WIDTHS = logspace(-3, 1, 9); % the kernel widths
   8 N = 100; % number of data points
   9 
  10 % set up some empty vectors to collect the information
  11 BIAS2 = zeros(length(WIDTHS), 1);
  12 VARIANCE = zeros(length(WIDTHS), 1);
  13 ERR = zeros(length(WIDTHS), 1);
  14 
  15 % plot individual means and variances for 9
  16 % different widths
  17 figure(1)
  18 for WI = 1:length(WIDTHS)
  19   W = WIDTHS(WI)
  20   
  21   % collect predictions for ITERS many resamples
  22   Y = zeros(M, ITERS);
  23   for I = 1:ITERS
  24     [DX, DY] = generate_data(N);
  25     Y(:, I) = kernel_regression(W, DX, DY, X);
  26   end
  27 
  28   % compute mean and variance
  29   ME = mean(Y, 2);
  30   VA = variance_per_point(Y);
  31 
  32   % plot the result
  33   subplot(3, 3, WI)
  34   plot(DX, DY, '.', ... % the last data set
  35        X, f(X), 'b-', ... % the true function 
  36        X, ME, 'r-', ... % the mean of the predictions
  37        X, ME + 2*sqrt(VA), 'r--', ... % +- 2*sqrt(variance)
  38        X, ME - 2*sqrt(VA), 'r--', ...
  39        X, Y(:, end), 'g-'); % the last learned function
  40 
  41   % compute mean squared bias, variance, and error
  42   BIAS2(WI) = mean(squared_bias_per_point(Y, f(X)));
  43   VARIANCE(WI) = mean(VA);
  44   ERR(WI) = mean(squared_error_per_point(Y, f(X)));
  45   title(sprintf('width = %f\n bias^2 = %f, var = %f', ...
  46                 W, BIAS2(WI), VARIANCE(WI)));
  47   drawnow
  48 end
  49 
  50 % plot squared bias, variance, and error
  51 figure(2)
  52 semilogx(WIDTHS, BIAS2, WIDTHS, VARIANCE, WIDTHS, ERR)
  53 legend('bias^2', 'variance', 'error')
  54 
  55 % Generate data 
  56 function [X, Y] = generate_data(N)
  57 X = sort(rand(N, 1) * 8 - 4);
  58 Y = f(X) + 0.5 * randn(N, 1);
  59 
  60 % the true function (sinc function + linear offset)
  61 function Y = f(X)
  62 Y = sin(4*X)./(4*X) + X/2;
  63 
  64 % Compute kernel regression estimates
  65 function YE = kernel_regression(W, X, Y, XE)
  66 K = rbfkern(W, X, XE);
  67 YE = Y'*K ./ sum(K);
  68 
  69 % Compute a Gaussian kernel
  70 function K = rbfkern(W, X1, X2)
  71 D = pwdist(X1, X2);
  72 K = exp(-D / W);
  73 
  74 % Compute all pairwise distances
  75 function D = pwdist(X, Y)
  76 D = size(X, 2);
  77 N = size(X, 1);
  78 M = size(Y, 1);
  79   
  80 XX = sum(X.*X, 2);
  81 YY = sum(Y.*Y, 2);
  82 D = repmat(XX, 1, M) + repmat(YY', N, 1) - 2*X*Y';
  83 
  84 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  85 % Your solutions below
  86 
  87 % 1. Compute the variance per point
  88 %
  89 % Arguments:
  90 %  Y: matrix where the columns are the predictions for a single
  91 %     version of the data set
  92 function V = variance_per_point(Y)
  93 % ...
  94 
  95 % 2. Compute the squared bias per point
  96 %
  97 % Arguments:
  98 %  Y: see 1
  99 %  F: a column vector of the true function
 100 function B2 = squared_bias_per_point(Y, F)
 101 % ...
 102 
 103 % 3. Compute the squared error per point.
 104 %
 105 % Arguments:
 106 %  Y: see 1
 107 %  F: see 2
 108 function E = squared_error_per_point(Y, F)
 109 % ...