lambdas=exp(-3:.5:2);
sigma=sqrt(1/10);
[xTest,yTest,yTrue]=genData(1000,sigma);

M=10; % model order for regression
XTestM=degexpandScale(xTest,M,ones(1,M+1));

K=200; % number of random training sets to try
yHat=zeros(K,1000);

clear testError biasSquared variance;

N = 20; % size of the random training sets

for l=1:length(lambdas)
  lambda=lambdas(l);

  for k=1:K
    % generate a random training set
    [x,y]=genData(N,sigma);
    X=degexpandScale(x,M,ones(1,M+1));
    
    % ridge regression with parameter lambda
    wHat =  ???? % fill this in
    yHat(k,:)=(XTestM*wHat)';
  end

  % done with this lambda; compute bias/variance/test error
  
  % average (over random training sets) prediction for each x
  yAvg=mean(yHat);

  % test error: averaged over test points AND over random training sets
  testError(l) = ???? % fill this in
  
  % squared bias: averaged over test points
  biasSquared(l) = ???? % fill this in
  
  % variance: averaged over test points
  variance(l) = ???? % fill this in
  
  fprintf(2,'Done with lambda=%.2f\n',lambda);
end

figure(1);clf;
semilogx(lambdas,biasSquared,'b');
hold on;
semilogx(lambdas,variance,'r');
semilogx(lambdas,biasSquared+variance,'m--');
semilogx(lambdas,testError,'k');
legend('bias^2','variance','bias^2+variance','test error');