Michael Chuah's page for MAS.864 The Nature of Mathematical Modelling

mcx "at" mit "dot" edu

Home

Function Fitting

Scanned work

MATLAB code output

11-1
 
errA =
 
    0.0419
 
 
errB =
 
    0.1145
 
 
errAb =
 
    0.0419
 
 
errBb =
 
    0.1150
 
 
errAi =
 
    0.0400
 
 
errBi =
 
    0.1207

MATLAB code used

%% Meng Yee CHUAH
% Week 8 - Function Fitting
close all; clear all; clc;
 
%% 11.1 - SVD
% Part (a)
x = rand(100,1);
g = randn(100,1)*0.5;
y = 2+3*x+g;
 
A = [ones(100,1) x];
[U W V] = svd(A,0);
a = V/W*U'*y;
 
xx = 0:0.1:1;
yy = a(1)+a(2)*xx;
figure('Units','inches','Position',[0 0 6 4.5],'PaperPosition',[0 0 6 4.5]); 
plot(x,y,'.',xx,yy); title('SVD of y=2+3x+\xi');
saveas(gcf,'11-1.png');
saveas(gcf,'11-1.fig');
 
errA = ((V(1,1)^2)/(W(1,1)^2))+((V(1,2)^2)/(W(2,2)^2))
errB = ((V(2,1)^2)/(W(1,1)^2))+((V(2,2)^2)/(W(2,2)^2))
 
% Part (b) - Bootstrap sampling to generate 100 data sets
errAb = zeros(100,1);
errBb = zeros(100,1);
for j = 1:100
    yb = randsample(y,100,true);
    xb = zeros(100,1);
    for i = 1:100
        xb(i) = x(find(y==yb(i),1));
    end
 
    Ab = [ones(100,1) xb];
    [Ub Wb Vb] = svd(Ab,0);
    ab = Vb/Wb*Ub'*yb;
 
%     xxb = 0:0.1:1;
%     yyb = ab(1)+ab(2)*xxb;
%     figure('Units','inches','Position',[0 0 6 4.5],'PaperPosition',[0 0 6 4.5]); 
%     plot(xb,yb,'.',xxb,yyb); title('Bootstrap');
 
    errAb(j) = ((Vb(1,1)^2)/(Wb(1,1)^2))+((Vb(1,2)^2)/(Wb(2,2)^2));
    errBb(j) = ((Vb(2,1)^2)/(Wb(1,1)^2))+((Vb(2,2)^2)/(Wb(2,2)^2));
end
mean(errAb)
mean(errBb)
 
% Part (c) - Ensemble of 100 independent data sets
for i = 1:100
    xi = rand(100,1);
    gi = randn(100,1)*0.5;
    yi = 2+3*xi+gi;
 
    Ai = [ones(100,1) xi];
    [Ui Wi Vi] = svd(Ai,0);
    ai = Vi/Wi*Ui'*yi;
 
%     xxi = 0:0.1:1;
%     yyi = ai(1)+ai(2)*xxi;
%     figure('Units','inches','Position',[0 0 6 4.5],'PaperPosition',[0 0 6 4.5]); 
%     plot(xi,yi,'.',xxi,yyi); title('Ensemble');
 
    errAi(i) = ((Vi(1,1)^2)/(Wi(1,1)^2))+((Vi(1,2)^2)/(Wi(2,2)^2));
    errBi(i) = ((Vi(2,1)^2)/(Wi(1,1)^2))+((Vi(2,2)^2)/(Wi(2,2)^2));
end
mean(errAi)
mean(errBi)