function nma_Math504_HW1(nBins,nSims,seed) %by Nasser Abbasi, HW1, Math 504 %calculate expected value of y by simulation %clear all; close all; % % C O N S T A N T S and P A R A M E T E R S % %nSims = 20000; %nBins = 50; %seed = 01010101; skip = round(0.01*nSims); a = 0; b = 1; str = 'HW1 Mathematics 504 CSUF spring 2008'; str =[str '\nsimulation[%d], std[%4.3f], mean[%6.5f], true mean[%5.4f]']; barWidth = 2/nBins; % % I N I T I A L I Z A T I O N % rand('twister',seed); box = zeros(1,2); y = zeros(nSims,1); figure(1); % % L O G I C % for i=1:nSims box(1)=rand; box(2)=2*box(1); boxSelected=((a + (b-a).*rand)<=.5)+1; y(i)=box(boxSelected); if y(i)<=1 %switch box if needed if boxSelected==1 boxSelected=2; else boxSelected=1; end y(i)=box(boxSelected); end [n,x]=hist(y(1:i),nBins); currentArea = barWidth*sum(n); if mod(i,skip)==0 bar(x,n/currentArea,'y'); %relative frequency title(sprintf(str,i,std(y),mean(y),15/16)); xlabel('x=final observed value'); ylabel('P(X=x)'); xtrue=[0 .5 .5 1 1 2]; ytrue=[.75 .75 .25 .25 .5 .5]; line(xtrue,ytrue,'Color','r','LineWidth',4); ylim([0,1]); legend('simulation PDF','exact PDF'); drawnow; pause(0.01) end end end
function nma_Math504_HW1_part2() %by Nasser Abbasi, HW1, part 2, Math 504 %calculate expected value of y by simulation %This scripts simulates the pdf of the observed value from %the following expeirment: % %pick a random number x from uniform[0,1], put this %number in a box, and put twice the number in a second %box. Next, pick one of these boxes by random, and look %at the number inside. Call this y. If the number is smaller than %one, then switch the box. Find the pdf of final y observed y, and %find the estimate of the mean of y. %Notice that 0<= y <= 2. clear all; close all; seed = 01010101; rand('twister',seed); % First do an initial estimate using simulation to estimate % population mean and standard deviation, and then use these to % obtain the needed sample size for the error level required [s,xBar] = initialEstimate(); err = 0.01 * xBar; sampleSize = ((1.96*s)/err)^2; % % C O N S T A N T S and P A R A M E T E R S % nBins = 50; skip = round(0.01*sampleSize); barWidth = 2/nBins; % % I N I T I A L I Z A T I O N % y = zeros(sampeSize,1); figure(1); set(0,'DefaultTextinterpreter','none'); h=title(''); axpos = get(gca,'pos'); extent = get(h,'extent'); set(gca,'pos',[axpos(1) axpos(2) axpos(3) axpos(4)-.45*extent(4)]); set(h,'VerticalAlignment','Middle'); % % L O G I C % for i=1:sampleSize y(i) = makeAnObservation(); % currentArea = barWidth*sum(n); if mod(i,skip)==0 generateOneFrame(y(1:i),nBins,sampleSize); end end generateFinalResult(y,nBins); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [s,xbar]=initialEstimate() sampleSize = 20000; y = zeros(sampleSize,1); for i = 1:sampleSize y(i) = makeAnObservation() end s = std(y); xbar = mean(y); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y = makeAnObservation() box = zeros(1,2); box(1) = rand; box(2) = 2*box(1); boxSelected = (rand<.5)+1; %pick a box by random y =box(boxSelected); if y <=1 %switch box if needed if boxSelected == 1 boxSelected = 2; else boxSelected = 1; end y = box(boxSelected); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function generateOneFrame(y,nBins,sampleSize) currentFrameNumber = length(y); firstLineTitle = 'HW1 Mathematics 504, part(2) CSUF spring 2008'; secondLineTitle = 'simulation[$%d$], $s=%6.5f$, $\\bar{X}=%9.8f$, $\\mu=%5.4f$'; [n,x]=hist(y,nBins); bar(x,n/sampleSize,'y'); %relative frequency ls=sprintf(secondLineTitle,currentFrameNumber,std(y),mean(y),15/16); h=title(char(firstLineTitle,ls),'fontsize',12,'interpreter','latex'); xlabel('x=final observed value'); ylabel('P(X=x)'); xtrue=[0 .5 .5 1 1 2]; ytrue=[.75 .75 .25 .25 .5 .5]; line(xtrue,ytrue,'Color','r','LineWidth',2,'LineStyle','--'); ylim([0,1]); legend('simulation PDF','exact PDF'); drawnow; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function generateFinalResult(y,nBins) sampleSize = length(y); firstLineTitle = 'HW1 Mathematics 504, CSUF spring 2008'; secondLineTitle = 'simulation[$%d$], $s=%6.5f$, $\\bar{X}=%9.8f$, $\\mu=%5.4f$'; % % Compare estimated mean with theortical mean % estimatedMean = mean(y); standardError = 1.96*std(y)/sqrt(sampleSize); lss=sprintf('95 perecent confidence interval for mean is (%6.5f ... %6.5f)',... estimatedMean-standardError,estimatedMean+standardError); ls=sprintf(secondLineTitle,sampleSize,std(y),mean(y),15/16); title({firstLineTitle,ls,lss},'fontsize',12,'interpreter','latex'); [n,x]=hist(y,nBins); bar(x,n/sampleSize,'y'); %relative frequency ls=sprintf(secondLineTitle,sampleSize,std(y),mean(y),15/16); h=title(char(firstLineTitle,ls),'fontsize',12,'interpreter','latex'); xlabel('x=final observed value'); ylabel('P(X=x)'); xtrue=[0 .5 .5 1 1 2]; ytrue=[.75 .75 .25 .25 .5 .5]; line(xtrue,ytrue,'Color','r','LineWidth',2,'LineStyle','--'); ylim([0,1]); legend('simulation PDF','exact PDF'); drawnow; end
% file nma_Math504_HW1_part1_as_script %by Nasser Abbasi, HW1, Math 504 %This scripts simulates the pdf of the observed value from %the following expeirment: % %pick a random number x from uniform[0,1], put this %number in a box, and put twice the number in a second %box. Next, pick one of these boxes by random, and look %at the number inside. Call this y. Find the pdf of %y. Notice that 0<= y <= 2. clear all; close all; % % C O N S T A N T S and P A R A M E T E R S % nSims = 20000; nBins = 50; seed = 01010101; %my seed to reproduce same results. skip = round(0.01*nSims); str = 'HW1 Mathematics 504, part 1. CSUF spring 2008'; str =[str '\nsimulation[%d], std[%4.3f], mean[%6.5f], true mean[%5.4f]']; barWidth = 2/nBins; % % I N I T I A L I Z A T I O N % rand('twister',seed); box = zeros(1,2); y = zeros(nSims,1); figure(1); % % L O G I C % for i=1:nSims box(1) = rand; box(2) = 2*box(1); boxSelected = (rand<.5)+1; y(i)=box(boxSelected); [n,x]=hist(y(1:i),nBins); currentArea = barWidth*sum(n); if mod(i,skip)==0 bar(x,n/currentArea,'y'); %relative frequency title(sprintf(str,i,std(y),mean(y),3/4)); xlabel('x=observed value'); ylabel('P(X=x)'); xtrue=[0 1 1 2]; ytrue=[.75 .75 .25 .25]; line(xtrue,ytrue,'Color','r','LineWidth',2,'LineStyle','--'); ylim([0,1]); legend('simulation PDF','exact PDF'); drawnow; %pause(0.01) end end title(sprintf(str,nSims,std(y),mean(y),3/4));
function nma_Math504_HW1_part2() %function nma_Math504_HW1_part2() % %This function simulates the pdf of the observed value from %the following expeirment: % %pick a random number x from uniform[0,1], put this %number in a box, and put twice the number in a second %box. Next, pick one of these boxes by random, and look %at the number inside. Call this y. If the number is smaller than %one, then switch the box. Find the pdf of final y observed y, and %find the estimate of the mean of y. %Notice that 0<= y <= 2. %by Nasser Abbasi, HW1, part 2, Math 504 clear all; close all; seed = 01010101; rand('twister',seed); % First do an initial estimate using simulation to estimate % population mean and standard deviation, and then use these to % obtain the needed sample size for the error level required [s,xBar] = initialEstimate(); err = 0.01 * xBar; sampleSize = round(((1.96*s)/err)^2); fprintf('s=%f, xBar=%f\n',s,xBar); % % C O N S T A N T S and P A R A M E T E R S % nBins = 50; skip = round(0.01*sampleSize); %for simulation, skip frames % % I N I T I A L I Z A T I O N % y = zeros(sampleSize,1); figure(1); set(0,'DefaultTextinterpreter','none'); h=title({'','',''}); axpos = get(gca,'pos'); extent = get(h,'extent'); set(gca,'pos',[axpos(1) axpos(2) axpos(3) axpos(4)-.20*extent(4)]); set(h,'VerticalAlignment','Middle'); % % L O G I C % for i=1:sampleSize y(i) = makeAnObservation(); if mod(i,skip)==0 generateOneFrame(y(1:i),nBins); end end generateFinalResult(y,nBins); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [s,xbar]=initialEstimate() sampleSize = 20000; y = zeros(sampleSize,1); for i = 1:sampleSize y(i) = makeAnObservation(); end s = std(y); xbar = mean(y); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y = makeAnObservation() box = zeros(1,2); box(1) = rand; box(2) = 2*box(1); boxSelected = (rand<.5)+1; %pick a box by random y = box(boxSelected); if y <=1 %switch box if needed if boxSelected == 1 boxSelected = 2; else boxSelected = 1; end y = box(boxSelected); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function generateOneFrame(y,nBins) barWidth = 2/nBins; sampleSize = length(y); firstLineTitle = 'HW1 Mathematics 504, part(2) CSUF spring 2008'; secondLineTitle = 'simulation[$%d$], $s=%6.5f$, $\\bar{X}=%9.8f$, $\\mu=%5.4f$'; [n,x] = hist(y,nBins); currentArea = barWidth*sum(n); bar(x,n/currentArea,'y'); %relative frequency ls=sprintf(secondLineTitle,sampleSize,std(y),mean(y),15/16); title(char(firstLineTitle,ls),'fontsize',12,'interpreter','latex'); xlabel('x=final observed value'); ylabel('P(X=x)'); xtrue = [0 .5 .5 1 1 2]; ytrue=[.75 .75 .25 .25 .5 .5]; line(xtrue,ytrue,'Color','r','LineWidth',2,'LineStyle','--'); ylim([0,1]); legend('simulation PDF','exact PDF'); drawnow; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function generateFinalResult(y,nBins) barWidth = 2/nBins; sampleSize = length(y); firstLineTitle = 'HW1 Mathematics 504, CSUF spring 2008'; secondLineTitle = 'simulation[$%d$], $s=%6.5f$, $\\bar{X}=%9.8f$, $\\mu=%5.4f$'; % % Compare estimated mean with theortical mean % estimatedMean = mean(y); standardError = 1.96*std(y)/sqrt(sampleSize); lss=sprintf('95 perecent confidence interval for mean is (%6.5f ... %6.5f)',... estimatedMean-standardError,estimatedMean+standardError); ls=sprintf(secondLineTitle,sampleSize,std(y),mean(y),15/16); [n,x]=hist(y,nBins); currentArea = barWidth*sum(n); bar(x,n/currentArea,'y'); %relative frequency title({firstLineTitle,ls,lss},'fontsize',12,'interpreter','latex'); xlabel('x=final observed value'); ylabel('P(X=x)'); xtrue=[0 .5 .5 1 1 2]; ytrue=[.75 .75 .25 .25 .5 .5]; line(xtrue,ytrue,'Color','r','LineWidth',2,'LineStyle','--'); ylim([0,1]); legend('simulation PDF','exact PDF'); drawnow; end
function [no,xo] = rhist(varargin) %RHIST Relative Histogram. % N = HIST(Y) bins the elements of Y into 10 equally spaced containers % and returns the relative frequency of elements in each container. If Y is a % matrix, RHIST works down the columns. % % N = RHIST(Y,M), where M is a scalar, uses M bins. % % N = RHIST(Y,X), where X is a vector, returns the relative freqency of Y % among bins with centers specified by X. The first bin includes % data between -inf and the first center and the last bin % includes data between the last bin and inf. Note: Use HISTC if % it is more natural to specify bin edges instead. % % N = RHIST(Y,M,Any_Character) returns relative frequency density of % Y among bins.Any_Character is the any character inside single quotation % or any numeric value. % You can omit second optional argument using single quotation % i.e. N = RHIST(Y,'',Any_Character) returns relative frequency density % for 10 bins. % It is to be noted that sum(N)equals unity for relative frequency % while area under curve for relative frequency density equals unity. % Note that as size(Y,1) and M increases relative frequency density is % close to probability density for continous random variable. % % [N,X] = RHIST(...) also returns the position of the bin centers in X. % % RHIST(...) without output arguments produces a histogram of relative % frequency or relative frequency densisty bar plot of the results. % The bar edges on the first and last bins may extend to cover the min % and max of the data unless a matrix of data is supplied. % % RHIST(AX,...) plots into AX instead of GCA. % % Class support for inputs Y, X: % float: double, single % % See also HIST. % Copyright 2004-2005 Durga Lal Shrestha. % $Revision: 1.0.0 $ $Date: 2005/6/20 14:30:00 $ % Parse possible Axes input error(nargchk(1,inf,nargin)); [cax,args,nargs] = axescheck(varargin{:}); y = args{1}; if nargs == 1 x = 10; elseif nargs == 2 x = args{2}; else if isempty(args{2}) x = 10; else x = args{2}; end end [m,n] = size(y); [nn,x]=hist(y,x); % frequency nn = nn./m; % relative frequency % relative frequency density if nargs == 3 binwidth = x(2)-x(1); nn = nn./binwidth; end if nargout == 0 if ~isempty(cax) bar(cax,x,nn,[min(y(:)) max(y(:))],'hist'); else bar(x,nn,[min(y(:)) max(y(:))],'hist'); end xlabel('y') if nargs == 3 ylabel('relative frequency density') else ylabel('relative frequency') end else no = nn; xo = x; end
l.174 — TeX4ht warning — \SaveEverypar’s: 3 at \begindocument and 4 \enddocument —