% SAMPLEERROR - Computes the error for three scenarios: a single uniform
% sample with Gaussian noise, bootstrapping that sample a number of times, and
% a number of independent samples.
%
%          samples - number of samples to produce
%               sd - standard deviation of the Gaussian noise
%             sets - number of bootstraps or independent sets to compute
%
% Peter Torpey
%


function sampleError(samples, sd, sets)
    xx = rand(samples, 1);
    nn = sd .* randn(samples, 1);
    aa = [2 3];
    gv = sd ^ 2;
    yy = aa(1) + aa(2) .* xx + nn;
    
    disp('a. Single Error');
    A = [ones(samples, 1) xx];
    [U, W, V] = svd(A);
    cc = V * pinv(W) * U' * yy;
    w = diag(W);
    e = sum((V(:, 1) .^ 2) ./ (w .^ 2));
    disp(sprintf('\ta Error: %0.5g\n', e(1)));
    e = sqrt(e * gv);
    disp(sprintf('\ta Value range: %0.5g - %0.5g\n', cc(1) - e(1), cc(1) + e(1)));
    e = sum((V(:, 2) .^ 2) ./ (w .^ 2), 2);
    disp(sprintf('\tb Error: %0.5g\n', e(2)));
    e = sqrt(e * gv);
    disp(sprintf('\tb Value range: %0.5g - %0.5g\n', cc(2) - e(2), cc(2) + e(2)));
    
    disp('b. Bootstrap Error');
    e1 = 0;
    e2 = 0;
    cc1 = 0;
    cc2 = 0;
    for i = 1:sets
        ss = xx(ceil(samples .* rand(samples, 1)));
        yy = aa(1) + aa(2) .* ss + nn;
        A = [ones(samples, 1) ss];
        [U, W, V] = svd(A);
        cc = V * pinv(W) * U' * yy;
        w = diag(W);
        e1 = e1 + sqrt(sum((V(:, 1) .^ 2) ./ (w .^ 2)) * gv);
        e2 = e2 + sqrt(sum((V(:, 2) .^ 2) ./ (w .^ 2)) * gv);
        cc1 = cc1 + cc(1);
        cc2 = cc2 + cc(2);
    end
    e1 = e1 / samples;
    e2 = e2 / samples;
    cc1 = cc1 / samples;
    cc2 = cc2 / samples;
    disp(sprintf('\ta Error: %0.5g\n', e1));
    disp(sprintf('\ta Value range: %0.5g - %0.5g\n', cc1 - e1, cc1 + e1));
    disp(sprintf('\tb Error: %0.5g\n', e2));
    disp(sprintf('\tb Value range: %0.5g - %0.5g\n', cc2 - e2, cc2 + e2));
    
    disp('c. Independent Sets Error');
    e1 = 0;
    e2 = 0;
    cc1 = 0;
    cc2 = 0;
    for i = 1:sets
        xx = rand(samples, 1);
        yy = aa(1) + aa(2) .* xx + nn;
        A = [ones(samples, 1) xx];
        [U, W, V] = svd(A);
        cc = V * pinv(W) * U' * yy;
        w = diag(W);
        e1 = e1 + sqrt(sum((V(:, 1) .^ 2) ./ (w .^ 2)) * gv);
        e2 = e2 + sqrt(sum((V(:, 2) .^ 2) ./ (w .^ 2)) * gv);
        cc1 = cc1 + cc(1);
        cc2 = cc2 + cc(2);
    end
    e1 = e1 / samples;
    e2 = e2 / samples;
    cc1 = cc1 / samples;
    cc2 = cc2 / samples;
    disp(sprintf('\ta Error: %0.5g\n', e1));
    disp(sprintf('\ta Value range: %0.5g - %0.5g\n', cc1 - e1, cc1 + e1));
    disp(sprintf('\tb Error: %0.5g\n', e2));
    disp(sprintf('\tb Value range: %0.5g - %0.5g\n', cc2 - e2, cc2 + e2));