Matlab code is
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
To run, save it to your Matlab working directory and type the command nmaHW4math504() from
the console.
l.389 — TeX4ht warning — \SaveEverypar’s: 3 at \begindocument and 4 \enddocument —