function [] = ConfEllipsoid_toy
close all;
seed = 1;
rng(seed);
n = 2e4;

x_i = (0:0.1:10);
modelfunc = @(x) 2*x+2;

y_biased = modelfunc(x_i)+randn(1, numel(x_i));
y = modelfunc(x_i);
figure(1);
hold on;
plot(x_i, y_biased, '.');
plot(x_i, y, '-');

beta0 = [1,1];
modelf = @(beta,x) beta(1)*x+beta(2);
[beta_opt,~,~,Cov_beta_opt,~,~] = nlinfit(x_i, y_biased, modelf, beta0);
figure(1);
plot(x_i, modelf(beta_opt, x_i), '-k');

betas = zeros(n,numel(beta_opt));
invCov_beta_opt = inv(Cov_beta_opt);
for i = 1:n
    beta_current = mvnrnd(beta_opt, Cov_beta_opt);
    while (beta_current-beta_opt)*invCov_beta_opt*(beta_current-beta_opt)' > chi2inv(0.9, length(Cov_beta_opt))
        beta_current = mvnrnd(beta_opt, Cov_beta_opt);
    end
    betas(i,:) = beta_current;
end

figure(2);
hold on;
plot(betas(:,1), betas(:,2), '.');
plot(beta_opt(1), beta_opt(2), 'ok');
min_betas = min(betas);
max_betas = max(betas);
plot(min_betas(1):0.0001:max_betas(1), min_betas(2), '-k', 'LineWidth', 3);
plot(min_betas(1):0.0001:max_betas(1), max_betas(2), '-k', 'LineWidth', 3);
plot(min_betas(1), min_betas(2):0.0001:max_betas(2), '-k', 'LineWidth', 3);
plot(max_betas(1), min_betas(2):0.0001:max_betas(2), '-k', 'LineWidth', 3);
xlim([min_betas(1)-0.05*max_betas(1) max_betas(1)+0.05*max_betas(1)]);
ylim([min_betas(2)-0.05*max_betas(2) max_betas(2)+0.05*max_betas(2)]);




figure(3);
beta1 = (min_betas(1)-0.05*max_betas(1):0.01:max_betas(1)+0.05*max_betas(1));
beta2 = (min_betas(2)-0.05*max_betas(2):0.01:max_betas(2)+0.05*max_betas(2));
[BETA1,BETA2] = meshgrid(beta1,beta2);
Z = mvnpdf([BETA1(:) BETA2(:)],beta_opt,Cov_beta_opt);
Z = reshape(Z, length(beta2), length(beta1));
surf(beta1,beta2,Z, 'EdgeColor', 'none');
grid on;
hold on;

end