% Script file for Math 305 - lab assignment 6


%  Text 1.7: 2
% weights of female black bears in age classes 1 to 15 are

W=[35 45 55 68 70 71 75 79 82 81 80 78 99 99 82];
A=1:15;

%  Fit the data with a straight line and determine the goodness of fit
n=length(A)
a=sum(A); aa=sum(A.^2); w=sum(W); aw=sum(A.*W);
M=[n, a; a , aa]  % linear regression matrix
d=[w ; aw]               % data about weights
p=M\d           % solve for p; entries are y-intercept and slope
xx=0:.1:15;
yy=p(1)+p(2)*xx;      % these are the y-coords of the regression line
predW=p(1)+p(2)*A;  % these are the predicted weights at each age
plot(A,W,'.',xx,yy,'-',A,predW,'*','MarkerSize',15,'LineWidth',1.5)
%  The dots represent the data points, the solid line is the regression
%  line; the stars are the predicted values;  
%  now we calculate the goodness of fit.
SSE = sum((W-predW).^2)  % error sum of sqares
SST = sum((W-mean(W)).^2)    %  total sum of squares
rho2 = 1-SSE/SST
rho = sqrt(rho2)

%  Instead we can compute the correlation coefficient by:
corrcoef(A,W)
%  There is also a built in function polyfitx,y,n) for computing the 
%  coefficients of the best fitting polynomial of degree n, and 
%  a function polyval for evaluting the polynomial at points

p=polyfit(A,W,1)  % compare to what you found above
plot(A,W,'.',xx,polyval(p,xx),'-','MarkerSize',15,'LineWidth',1.5)

% Repeat the two lines replacing degree 1this with 
% degree 2, and then degree 3

%  Simulation
% Read about Matlab's random number generators by command:

help rand

x=1:100;
y=rand(1,100);
plot(x,y)

% Heads or Tails.  Flip a coin 40 times.  Alan wins one penny if heads
% loses one penny if tails.  X = Alan's profit; possible values are
% even numbers from -40 to 40.  Y=X/2 + 20 has values 0, 1, ... 40,
% and has the bin(40,.5) pmf.  So X=2(Y-20)

y=0:40; x=2*(y-20);
% next we give a very simple but clumsy construction of binomial pmf
bin=factorial(40)./(factorial(y).*factorial(40-y).*2^(40));
plot(x,bin)

% simulate the probability by using 100 repetiions of flipping a coin 40
% times

winningslist=[];leadslist=[];
for i=1:100
    winnings=0;  % don't forget these semicolons
    leads=0;
    for j=1:40
        lastwin=winnings;
        if rand()>.5
            winnings = winnings+1;
        else winnings = winnings -1;
        end
        if winnings >0 | (winnings == 0  & lastwin>0)
            leads=leads+1;
        end
    end
    winningslist=[winningslist,winnings];
    leadslist=[leadslist,leads];
end

hist(winningslist,x)  % histogram of winnings
% to compare to the binomial graph we have to scale up
% the vertical axis of the pmf 

hold on
 plot(x,100*bin)
hold off
 % repeat the simulation using 1000 trials instead of 100;
 % change the 100 in the for loop and in the plot to 1000
 
 % now make a histogram of the leadslist
 %  the leads can be from 0 to 40
 %  which values are least likely and which most?