clear all
close all
format compact
clc

Name       = 'Alex Bielajew'; % Substitute your name
Assignment = 3;


Problem = input('Which problem number [0 1 2 3 4 ] would you like to example? ');

if ((floor(Problem) ~= Problem) || (Problem < 0) || (Problem > 4))
   disp('Input must be integer an between 0 and 4 incusive. Try again.')
   return
end

N = 10000 % Number of samples. Change this to see effect on statistics

sprintf('This is %s''s solution for Problem %g.%g:',Name, Assignment, Problem);

if (Problem == 1)

elseif (Problem == 2)

elseif (Problem == 3)

elseif (Problem == 4)

elseif (Problem == 0) % p(x) = 2*x for 0 <= x <= 1

   disp('This is the example problem')
   
   M = 64;    % Mesh size of the plot

   figure1Handle = figure(1);
   set(figure1Handle,'Position',[222 493 568 410])
   hold on
   x  = [];
   x2 = [];
   for i = 1:N
      r        = rand;
      x(i)     = sqrt(r);
      x2(i)    = r;
      xbar(i)  = sum(x)/i;
      x2bar(i) = sum(x2)/i;
      s2x(i)   = sum((x - xbar).^2)/(i-1);
   end
   plot1Handle = plot( ...
      log10(1:N),xbar,'k.',        ...
      log10([1,N]),[2/3,2/3],'k-', ...
      log10(1:N),x2bar,'b.',       ...
      log10([1,N]),[1/2,1/2],'b-', ...
      log10(1:N),s2x,'r.',         ...
      log10([1,N]),[1/18,1/18],'r-'...
   );
   legendHandle = legend(                               ...
      '$\langle x\rangle$  ' ,'theory: $\mu  = \int_0^1{\rm d}x\,x\, p(x) = \frac{2}{3}$ ', ...
      '$\langle x^2\rangle$' ,'theory: $\int_0^1{\rm d}x\,x^2\, p(x) = \frac{1}{2}$ ', ...
      '$\sigma^2_x$        ' ,'theory: $\sigma^2_x = \int_0^1{\rm d}x\,(x - \mu)^2\, p(x) = \frac{1}{18}$'  ...
  );
   set(legendHandle,'Interpreter','LaTex','FontSize',12,'Position', [0.3945 0.2182 0.5307 0.3217])
   xlabel('$\log_{10}(N)$','Interpreter','LaTex','FontSize',20)
   ylabel('$\langle x\rangle$, $\langle x^2\rangle$, and $\sigma^2_x$',...
          'Interpreter','LaTex','FontSize',20)
   title('{\bf $p(x) = 2 x;~0\le x\le 1$: Convergence of statistics}','Interpreter','LaTex','FontSize',18)
   % Sort them into bins
   xbar   = sum(x)/N
   x2bar  = sum(x2)/N
   s2x    = sum((x - xbar).^2)/(N-1)
   s2xbar = s2x/N
   
   figure(2) % This figure shows how faithfully the pdf is reproduced
   p = zeros(1,M); % Sizes the p vector, a row vector
   for (i = 1:N)
      index = ceil(M*x(i)); % indexes measurement datum into the display mesh
      p(index) = p(index) + 1;
   end
   praw = p; % Save for later for individual error bars
   p = M*p/N; % Normalization /N to normalize/history
             %               *M due to intrinsic (1/M) mesh size (i.e. dx)
   x = linspace(0,1,M); % Resizes the x vector for display   
   plot(x,p,'ko', [0,1],[0,2],'k')
   xlabel('$x$','Interpreter','LaTex','FontSize',20')
   ylabel('$p(x)$','Interpreter','LaTex','FontSize',20)
   title('{\bf Sampled and theoretical $p(x) = 2 x$}','Interpreter','LaTex','FontSize',18)
   legendHandle = legend('Monte Carlo','Theory');
   set(legendHandle,'Interpreter','LaTex','FontSize',12,'Position', [0.5803 0.2458 0.2272 0.1046])
   
   % Plot the individual error bars
   hold on % Keep the current plot
   hi = (M/N)*(praw + sqrt(praw));
   lo = (M/N)*(praw - sqrt(praw));
   for i = 1:M
      plot([x(i),x(i)],[lo(i),hi(i)],'k-') % Plots the error bars using Poisson statistics
   end
   
end