function nmaHW4math504()
%function nmaHW4math504()
%Solve problem #2 in second handout, Math 504
%spring 2008 CSUF
%by Nasser Abbasi  feb 8,2008

%Developed on MATLAB Version 7.4.0.287 (R2007a)
%Running on Win XP. Uses statistics toolbox

%
%  M A I N    C O N F I G U R A T I O N     S E C T I O N
%
params.D = 3;          % Diffusion parameter
params.T = 2;          % total running time
params.beta = 2;       % drift parameter

%
%
% These parameters are derived automatically from the above
%
params.trueVar  = 2*params.D*params.T;
params.trueStd  = sqrt(params.trueVar);
params.trueMean = params.beta*params.T;

%
%  internal program   C O N F I G U R A T I O N
%  These are configutation parameters for displying
%  and for setting sample size and number of steps
%  Adjust as needed
%

config.seed  = 12345;
config.nBins = 10;
config.sample_size = 5000;
config.max_number_of_steps = 20;

%
%determine the starting number of steps such that 'p' comes out to be
%less than 1. see report for derivation of this formula
%
config.starting_step = round(params.T*params.beta^2/(2*params.D)) +1;

%add the above number of steps to the starting step to
%obtain max number of steps

config.max_number_of_steps = config.starting_step+config.max_number_of_steps;

%
% set the number of steps to skip at each simulation else this
% will take too long to run
%
config.n = config.starting_step:5:config.max_number_of_steps;

%
% Rest are internal data structures to keep track of
% simulation data during runs
%
config.varianceVector = zeros(length(config.n),1);
config.rmsVector      = zeros(length(config.n),1);
config.delT           = zeros(length(config.n),1);
config.standarize     = 1;  %set this to zero if you do not standarize

%
% Determine the theoretical quantiles for the standard normal
% distribution to use in the plots generated
%
config.qp=norminv( ((1:config.sample_size)-.5)/config.sample_size,0,1);

%
% I N I T I A L I Z A T I O N
% Create the figure and seed the random number generator
%

makeFigure();
hold on;
rand('state',config.seed);

%
% Here we go, let start the fun part
%

for i = 1:length(config.n)
    config=simulate_one_walk( params, config, config.n(i), i );
end

end

%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%
function config=simulate_one_walk( params, config, number_of_steps, ...
    current_experiement_number )

%
% generate delt(T) and del(X) and p and q from the input
%
config.delT(current_experiement_number) = params.T/number_of_steps;
delX = sqrt(2*params.D*config.delT(current_experiement_number));
p    = 1/2 * (1 + params.beta*delX/(2*params.D));
q    = 1-p;

[normalizedPosition,position] = generate_distribution(p,q,delX,...
    number_of_steps,params,config);

config.varianceVector(current_experiement_number) = var(position);

if config.standarize
    pos=normalizedPosition;
else
    pos=position;
end

%
% Now that we have a distribution generated, lets find the rms error
%
[config,truePDF,estimatedPDF,estimatedFit,xForSimulation,xForNormal]=...
    getRMSerrorInCurrentPDF(pos,config,params,p,q,...
    current_experiement_number);

%
%Ok, we have all the data, lets make a plot
%
updatePlots(number_of_steps,normalizedPosition,...
    delX,params,config,truePDF,estimatedPDF,p,q,...
    xForNormal,current_experiement_number);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  This function generates a sample
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [normalizedPosition,position]=generate_distribution(p,q,...
    delX,number_of_steps,params,config)

%
% generate 2 arrays to hold the final positions in. One
% for standardized position and for not standarized position
% (was not sure which to use at one point, so I keep both)
%
position = zeros(config.sample_size,1);
normalizedPosition = position;

for i = 1:config.sample_size
    y = makeOneRandomWalk(p,q,number_of_steps);
    position(i) = sum(y)*delX;
    normalizedPosition(i) = (position(i)-params.trueMean) / params.trueStd;
end

end
%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%
function y=makeOneRandomWalk(p,q,number_of_steps)

y       = rand(number_of_steps,1);
y(y<=q) = -1;
y(y>q)  = 1;

end
%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%
function updatePlots(current_number_of_steps,...
    normalizedPosition,delX,params,config,truePDF,...
    estimatedPDF,p,q,xForNormal,current_experiement_number)

subplot(3,2,1:2);
str1='first step[$%d$] current step[$%d$] sampleSize[$%d$] $p=[%4.3f]$ $q=[%4.3f$] $\\beta t=[%4.3f]$ $D = [%d]$ $t = [%d]$';
cla;
bar(config.xout,estimatedPDF,'y');  %relative frequency

line1=sprintf(str1,config.starting_step,current_number_of_steps,...
    config.sample_size,p,q,params.trueMean,params.D,params.T);

title(char(line1),'fontsize',10,'interpreter','latex');
xlabel(sprintf('y'),'fontsize',10,'interpreter','latex');
ylabel('pdf(position)');
set(gca,'FontSize',8);
hold on;
plot(xForNormal,truePDF,'--r','LineWidth',2);

ymax = max(truePDF);
ylim([0,ymax+.5*ymax]);
if config.standarize
    xlim([-4,4]);
else
    xlim([params.trueMean-4*params.trueStd,...
        params.trueMean+4*params.trueStd]);
end
legend('current pdf','limit pdf');
drawnow;

%
% relative error in variance plot
%
subplot(3,2,3);
cla;
line([config.n(1) config.n(end)],[0  0]);
z=repmat(params.trueVar,current_experiement_number,1);
relativeErrorInVar = ...
    ((abs(z-config.varianceVector(1:current_experiement_number)))./z)*100;

line(config.n(1:current_experiement_number),relativeErrorInVar);
xlim([config.n(1) config.n(end)]);
ylim([0 80]);

line1 = sprintf('relative error in variance ($2DT$)');
line2 = 'true var$=%4.1f$, current var$=%4.1f$';
line2 = sprintf(line2,params.trueVar,config.varianceVector(current_experiement_number));
title(char(line1,line2),'fontsize',10,'interpreter','latex');
ylabel('relative error percentage','fontsize',10);
xlabel('number of steps (n)','fontsize',10);
%set(gca,'XTick',1:config.max_number_of_steps);
set(gca,'FontSize',8);
drawnow;

%
% RMS error plot
%
subplot(3,2,4);
cla;

plot(log2(config.delT(1:current_experiement_number)),...
    log2(config.rmsError(1:current_experiement_number)),'r-');
line1=sprintf('\\ \\ \\ \\ \\ \\  RMS error in PDF as $\\bigtriangleup{T}\\rightarrow{0}$');
line2=sprintf('$\\bigtriangleup{T}=%6.4f \\: \\bigtriangleup{X}=%6.4f \\:\\: rms=%6.5f$',...
    config.delT(current_experiement_number),delX, ...
    config.rmsError(current_experiement_number));
title(char(line1,line2),'fontsize',10,'interpreter','latex');

ylabel(sprintf('$log_2$(rms error)'),'fontsize',10,'interpreter','latex');
xlabel(sprintf...
    ('$log_2(\\bigtriangleup{T})$'),'fontsize',10,'interpreter','latex');

xlim([-6,1]);
ylim([-10,1]);
set(gca,'FontSize',8);
drawnow;

%
% quantile plot
%
subplot(3,2,5:6);
cla;
sp=sort(normalizedPosition);
plot(config.qp,sp,'r.');

xlabel('quantile of normal distribution');
ylabel('quantile of final position');
title('Quantile-Quantile plot','fontsize',10);
hold on;
line([-4 4],[-4 4]);
legend('data','normal distribution','Location','NorthWest');
set(gca,'FontSize',8);
xlim([-4,4]);
ylim([-4,4]);
drawnow;

end
%%%%%%%%%%%%%%%%%%
% Called to obtain the RMS error between the normal distibution
% and the current sample distribution
%
%%%%%%%%%%%%%%%%%%%%
function [config,truePDF,estimatedPDF,estimatedFit,...
    xForSimulation,xForNormal]=getRMSerrorInCurrentPDF(position,config,...
    params,p,q,current_experiement_number)

mu=mean(position);
stdd=std(position);

xForSimulation = linspace(mu-4*stdd,...
    mu+4*stdd,config.nBins);
[estimatedPDF,config.xout] = hist(position,xForSimulation);
config.binWidth = abs(  abs(config.xout(2))-abs(config.xout(1)));
currentArea     = config.binWidth*sum(estimatedPDF);
estimatedPDF    = estimatedPDF/currentArea;

if config.standarize
    stdd = 1;
    mu   = 0;
else
    stdd = params.trueStd;
    mu   = params.trueMean;
end

xForNormal = linspace(mu-4*stdd, mu+4*stdd, config.nBins);
truePDF    = pdf('Normal',xForNormal,mu, stdd);

diffPDF = truePDF-estimatedPDF;

config.rmsError(current_experiement_number) = ...
    norm(diffPDF)/sqrt(length(diffPDF));

estimatedFit = pdf('Normal',xForSimulation,params.trueMean, ...
    stdd*sqrt(4*p*q));

end
%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%
function makeFigure()

figure;
set(gcf,'Position',[200 100 700 600]);
set(gcf,'Resize','off')
set(0,'DefaultTextinterpreter','none');
axpos = get(gca,'pos');
h = title({'',''});
extent = get(h,'extent');

% position is [left, bottom, width, height];
set(gca,'pos',[axpos(1) axpos(2) axpos(3) axpos(4)-.3*extent(4)]);

end
