function central_limit_theorem
%built by keith regner 9/2016 with R2014b (8.4 OSX)
clear all
close all

%% description
% This demo illustrates the Central Limit Theorem. Consider an experiment 
% where n measurements are drawn from a non-normal population. The 
% measurements from this representative experiment are visualized in a 
% normalized histogram as shown in the top plot. The population from which 
% these measurements were drawn (solid black line) and the mean and 
% standard deviation of this population distribution are also shown in the 
% top plot. This experiment is repeated N times and the mean of the n 
% measurements from each of the N experiments are visualized in a 
% normalized histogram shown in the bottom plot. These data are then fit 
% with a normal distribution (solid black line) and the mean and standard 
% deviation of the resulting fitted normal distribution are indicated. 
% According to the Central Limit Theorem, the measurement means from each 
% experiment should follow a normal distribution with a mean equal to the 
% mean of the population distribution and standard deviation equal to the 
% population standard deviation over ?n if n ? 30. A comparison of the 
% population mean and standard deviation (indicated on the top plot) and 
% the fitted mean and fitted standard deviation multiplied with ?n 
% (indicated on the bottom plot) shows this fact to be quite accurate if
% n ? 30. Change the number of experiments N and the number of measurements 
% in each experiment n with the slider bars.  

%% slider parameters
b.num_samples = 100;                %number of experiments
b.sample_size = 10;                 %number of measurements
b.x = 0:.01:15;                     %range of original distribution
b.BW = 20/b.sample_size;            %binwidth for histogram of sample means
num_slide_min = 2;                  %num samples slider min
num_slide_max = 1002;               %num samples slider max
size_slide_min = 2;                 %sample size slider min
size_slide_max = 1002;              %sample size slider max

%% figure and axis formatting
b.fig = figure('Units','inches','Position',[5 .25 8 9],...
    'Name','central_limit_theorem','NumberTitle','off');
b.ax1 = axes('Units','inches','Position',[1 5.7 6 3],'FontSize',12,...
    'LineWidth',1.5,'XLim', [min(b.x) max(b.x)],'YLim',[0 0.5],'Box','on');
hold on
b.ax1.XLabel.String = 'X [a.u.]';
b.ax1.YLabel.String = 'PDF [a.u.]';
b.ax2 = axes('Units','inches','Position',[1 2 6 3],'FontSize',12,...
    'LineWidth',1.5,'XLim', [min(b.x) max(b.x)],'Box','on');
hold on
b.ax2.XLabel.String = 'X [a.u.]';
b.ax2.YLabel.String = 'PDF [a.u.]';

%% solution
%generate distribution
%lognormal
mu = 1;
sigma = 1;
b.ME = exp(mu + ((sigma^2)/2));
b.SD = sqrt(exp(2*mu + (sigma^2))*(exp(sigma^2) - 1));
b.pop_dist = pdf('Lognormal',b.x,mu,sigma);
b.pd  = makedist('Lognormal',mu,sigma);

%triangular
% a = 1;
% b = 14;
% c  = 15;
% b.ME = (a + b + c)/3;
% b.SD = sqrt(((a^2) + (b^2) + (c^2) - (a*b) - (a*c) - (b*c))/18);
% b.pop_dist = pdf('Triangular',b.x,'a',a,'b',b,'c',c);
% b.pd  = makedist('Triangular','a',1,'b',14,'c',15);

%uniform
% lower = 1;
% upper = 5;
% b.ME = (lower + upper)/2;
% b.SD = sqrt(((upper - lower)^2)/12);
% b.pop_dist = pdf('Uniform',b.x,lower,upper);
% b.pd  = makedist('Uniform',lower,upper);

%gamma
% a = 1;
% b = 1;
% b.ME = a*b;
% b.SD = sqrt(a*(b^2));
% b.pop_dist = pdf('gamma',b.x,a,b);
% b.pd  = makedist('gamma',a,b);

%select from distribution randomly and take average of each sample
%make arbitrary population distribution
b.experiment = FCN(b.sample_size,b.num_samples,b.pd);
b.sample_means = mean(b.experiment);

%create and plot histogram of sample means
b.hist = histogram(b.ax2,b.sample_means,'BinWidth',b.BW,...
    'Normalization','pdf');

%fit sample mean distribution with normal distribution
b.pd_fit = fitdist(b.sample_means','Normal');
b.X_fit = b.pd_fit.mu;              %mean of sample means
b.s_fit = b.pd_fit.sigma;           %standard deviation of sample means
b.fit_dist = pdf('Normal',b.x,b.X_fit,b.s_fit);

%plot fitted and population distributions
b.h_fit = plot(b.ax2,b.x,b.fit_dist,'Color',[0 0 0],...
    'LineWidth',1.5);
b.h_pop = plot(b.ax1,b.x,b.pop_dist,'Color',[0 0 0],'LineWidth',1.5);
b.ex_data = histogram(b.ax1,b.experiment(:,1),'BinWidth',b.BW,...
    'Normalization','pdf');
b.leg2 = legend(b.ax2...
    ,{'histogram of sample means','fitted normal distribution'},...
    'Location','NorthEast','Box','off','FontSize',12);
b.leg1 = legend(b.ax1...
    ,{'population distribution','histogram of one experiment'},...
    'Location','NorthEast','Box','off','FontSize',12);

%% add user control to numer of samples and sample size
b.txtboxMEpop = annotation('textbox','Units','inches','Position',...
    [5 7.5 4 .3],'String',...
    ['pop. mean = ' num2str(b.ME)],'FontSize',12,...
    'LineStyle','none');
b.txtboxSDpop = annotation('textbox','Units','inches','Position',...
    [5 7.2 4 .3],'String',...
    ['pop. STD = ' num2str(b.SD)],'FontSize',12,...
    'LineStyle','none');
b.txtboxMEnormfit = annotation('textbox','Units','inches','Position',...
    [5 3.8 4 .3],'String',...
    ['fitted mean = ' num2str(b.X_fit)],'FontSize',12,...
    'LineStyle','none');
b.txtboxSDnormfit = annotation('textbox','Units','inches','Position',...
    [5 3.5 4 .3],'String',...
    ['fitted STD = ' num2str(b.s_fit)],'FontSize',12,...
    'LineStyle','none');
b.txtboxSDovern = annotation('textbox','Units','inches','Position',...
    [5 3.2 4 .3],'String',...
    ['fitted STD*n^1^/^2 = ' num2str(b.s_fit*sqrt(b.sample_size))],...
    'FontSize',12,'LineStyle','none');
b.txtbox1 = annotation('textbox','Units','inches','Position',...
    [2.8 .85 4 .3],'String',...
    ['number of experiments (N) = ' num2str(b.num_samples)],'FontSize',12,...
    'LineStyle','none');
b.txtbox2 = annotation('textbox','Units','inches','Position',...
    [2.7 .05 4 .3],'String',...
    ['number of measurements (n) = ' num2str(b.sample_size)],'FontSize',12,...
    'LineStyle','none');

b.sl1 = uicontrol('Parent',b.fig,'Style','slide','Units','inches',...
    'Position',[1.5 1.2 5 .2],'Min',num_slide_min,'Max',num_slide_max,...
    'Value',b.num_samples,'Tag','slider1','SliderStep',[1/200 1/20],...
    'UserData',struct('NUM',b.num_samples,'SIZE',b.sample_size),...
    'Callback',{@slider1_call,b},'BackgroundColor',b.fig.Color);
b.label1a = uicontrol('Parent',b.fig,'Units','inches','Style',...
    'text','Position',[1.25 1.1 .2 .3],'String',num2str(num_slide_min),...
    'FontSize',12,'BackgroundColor',b.fig.Color);
b.label2a = uicontrol('Parent',b.fig,'Units','inches','Style',...
    'text','Position',[6.55 1.1 .4 .3],'String',num2str(num_slide_max),...
    'FontSize',12,'BackgroundColor',b.fig.Color);

b.sl2 = uicontrol('Parent',b.fig,'Style','slide','Units','inches',...
    'Position',[1.5 .4 5 .2],'Min',size_slide_min,'Max',size_slide_max,...
    'Value',b.sample_size,'Tag','slider2','SliderStep',[1/200 1/20],...
    'UserData',struct('NUM',b.num_samples,'SIZE',b.sample_size),...
    'Callback',{@slider2_call,b},'BackgroundColor',b.fig.Color);
b.label1b = uicontrol('Parent',b.fig,'Units','inches','Style',...
    'text','Position',[1.25 .3 .2 .3],'String',num2str(size_slide_min),...
    'FontSize',12,'BackgroundColor',b.fig.Color);
b.label2b = uicontrol('Parent',b.fig,'Units','inches','Style',...
    'text','Position',[6.55 .3 .4 .3],'String',num2str(size_slide_max),...
    'FontSize',12,'BackgroundColor',b.fig.Color);

set(findall(b.fig,'-property','FontUnits'),'FontUnits','normalized')
set(findall(b.fig,'-property','Units'),'Units','normalized')

end

%% slider1 for num of samples
function [] = slider1_call(varargin)
[H,b] = varargin{[1,3]};
h = findobj('Tag','slider2');
data = h.UserData;
slider_value = H.Value;
experiment = FCN(data.SIZE,slider_value,b.pd);
sample_means = mean(experiment);

pd_fit = fitdist(sample_means','Normal');
X_fit = pd_fit.mu;
s_fit = pd_fit.sigma;
fit_dist = pdf('Normal',b.x,X_fit,s_fit);
set(b.h_fit,'ydata',fit_dist)
data = struct('NUM',slider_value,'SIZE',data.SIZE);

%update plot and indicators
b.hist.Data = sample_means;
b.ex_data.Data = experiment(:,1);
b.hist.BinWidth = 20/data.SIZE;
b.ex_data.BinWidth = 20/data.SIZE;
set(b.txtbox1,'String',...
    ['number of experiments (N) = ' num2str(slider_value)])
set(b.txtboxMEnormfit,'String',...
    ['fitted mean = ' num2str(X_fit)]);
set(b.txtboxSDnormfit,'String',...
    ['fitted STD = ' num2str(s_fit)]);
set(b.txtboxSDovern,'String',...
    ['fitted STD*n^1^/^2 = ' num2str(s_fit*sqrt(data.SIZE))]);
H.UserData = data;
end

%% slider2 for sample size
function [] = slider2_call(varargin)
[H,b] = varargin{[1,3]};
h = findobj('Tag','slider1');
data = h.UserData;
slider_value = H.Value;
experiment = FCN(slider_value,data.NUM,b.pd);
sample_means = mean(experiment);

pd_fit = fitdist(sample_means','Normal');
X_fit = pd_fit.mu;
s_fit = pd_fit.sigma;
fit_dist = pdf('Normal',b.x,X_fit,s_fit);
set(b.h_fit,'ydata',fit_dist)
data = struct('NUM',data.NUM,'SIZE',slider_value);

%update plot and indicators
b.hist.Data = sample_means;
b.ex_data.Data = experiment(:,1);
b.hist.BinWidth = 20/slider_value;
b.ex_data.BinWidth = 20/data.SIZE;
set(b.txtbox2,'String',...
    ['number of measurements (n) = ' num2str(slider_value)])
set(b.txtboxMEnormfit,'String',...
    ['fitted mean = ' num2str(X_fit)]);
set(b.txtboxSDnormfit,'String',...
    ['fitted STD = ' num2str(s_fit)]);
set(b.txtboxSDovern,'String',...
    ['fitted STD*n^1^/^2 = ' num2str(s_fit*sqrt(slider_value))]);
H.UserData = data;
end

%% select data
function [experiment] = FCN(sample_size,num_samples,pd)
experiment = random(pd,[sample_size,num_samples]);
end

