% Unit 4 (MATLAB Course @ UAM)
% (c) Alexander Knebe

%====================================================================
% fit a Guassion to the distribution of the return values of randn()
%====================================================================
clear all, close all, clc

% obtain the normalised distribution
N     = 100000;
Nbins = 50;
data  = randn(N,1);
[p,x] = hist(data,Nbins);
p     = normalise(p,x);

% filter out zero entries in the distribution
i = find(p>0);
p = p(i);
x = x(i);

plot(x,p,'.'), hold on

% generate matrix Xij = dfi/daj for model: f = a1 - a2 x^2
dfda1 = ones(1,length(x));
dfda2 = -x.^2;
X     = [dfda1; dfda2]';


% obtain best-fit parameters via solution to linear system
a = inv(transpose(X)*X)*transpose(X) * log(p)'

% check if the two sigma values agree with each other
sigma1 = 1/(exp(a(1))^2*2*pi)
sigma2 = 1/sqrt(2*a(2))
std(data)

% plot best-fit
plot(x,exp(a(1)).*exp(-a(2).*x.^2),'r')
xlim([-5 5])
ylim([0 0.5])
xlabel('x')
ylabel('p(x)')
legend('distribution of randn() values','best-fit Gaussian')



