% Phd Computer Review Session: Matlab

% Vectors and Matrices
% 3x1 Vector (Column)
c = [1;2;3];
% 1x3 Vector (Row)
r = [1 2 3];
% 1x1 Matrix (scaler)
s = 7;
% Matrix: row vectors seperated by columns
m = [1 4 7 10; 2 5 8 11; 3 6 9 12];


% Matrix Operations
a = size(m);
g = [4 3 1];
h = g';                                             % returns the transpose of g
i = eye(3,2);                                       % identity matrices
ii = eye(2);
A = [1 1 1; 1 2 3; 1 3 6];
IA = inv(A);                                        % returns the inverse of A
d = det(A);
a = ones(2,3);                                      % matrices for ones
b = zeros(2,3);
c = zeros(size(A));
x = x(~isnan(x));                                   % remove unwanted NaN's from variables

A = [2 3; 4 5];
B = [2 1; 5 8];
r = [3 1 4];
s = 7;
C = A + B;                                          % Matrix addition

w = r * c;                                          % Matrix multiplication
%w = c * r;                                          
w = r * s;
C = A^2;
x = [3 1 4];
y = [5 6 4];
g = immultiply(x,y);                                % Term-by-term multiplication



% Stats and Random Variables
m = mean(r);
med = median(r);
s = std(r);
v = var(r);                                         % Variance of the data vector r
c = var(X);                                         % Cov matrix
int = iqr(r);                                       % Inter-quartile range

% random 5 x 5 matrix with uniform randomly variables
mr = rand(5);
mrn = randn(10);

% Putting the two together
M = 5000; N = 500; % Let?s choose M and N reasonably large.
A = rand(M,N); b = rand(N,1);
C = zeros(5000,1);
for i=1:5000
    for j=1:500
        C(i) = C(i) + A(i,j) * b(j);
    end
end
% Why is this stupid? What's a faster way to get the same result?

% Loading data and Plots
%load testdata.dat;                                 %   read data into testdata matrix
                                                    %       make sure there is no labtes (add
                                                    %       later if necessary)
month = testdata(:,1);                              %   copy first column of data into month
precip = testdata(:,2);                             %   and second column into precip
plot(month,precip,'o');                             %   plot precip vs. month with circles
xlabel('month of the year');                        % add axis labels and plot title
ylabel('mean precipitation (inches)');
title('Mean monthly precipitation at Portland International Airport');

% More plots
t = 0:0.1:10; % say we take samples every .1 seconds from 0 to 10 seconds
y = sin(2*pi*t/5); % make y a sine wave with period 5 seconds
plot(y); % the simplest way to plot; notice that the horizontal axis
% is not 0 to 10 as we would like, but rather corresponds to
% the indices of y, by default
plot(t,y); % specify both the horizontal and vertical quantities to
% plot against one another
plot(t,y,'go--'); % change the appearance of the plot with an optional
% third argument, a string that can specify color
% (here, green), point appearance (circles), and
% line style (dashed); look in ?help plot? for more
% options
x = cos(2*pi*t/5);
plot(x,y); % parametric plots are equally possible
plot(t,x,'g-.',t,y,'r*'); % more than one plot can be put on the same
% axes, by adding more triplets of arguments
% to the plot command
plot(t,x,'k:'); % note that making a new plot by default clears the
% old one
hold on; % after this command, future plots will appear on the same
% axes without clearing them
plot(t,y,'m-.');
hold off; % and now future plots will clear the axes again
plot(x,y);
semilogx(x+2,y+2); % make the x-axis logarithmic
semilogy(x+2,y+2); % or the y-axis
loglog(x+2,y+2); % or both axes



% Linear Regression
N = length(precip);
X = [ones(N,1), month];
[bhat bint R Rint Stats] = regress(precip,X);
% Type help regress to find out what the returned data are.

% Non-parametric bootstrap: re-sampling with replacement from the errors
Yhat = precip-R;                                    %   Calculate predicted Y values
bootb = bootstrp(1000, @(x) regress(Yhat+x,X),R);   %   Bootstrap residuals R
bootb = sort(bootb);                                %   Rank the coefficients from low to high
bootCI = [bootb(25,:); bootb(975,:)];               %   Report 5% confidence intervals


% Regression with lags
lagY = [NaN; precip(1:end-1)];
X = [X lagY];
[bhat bint R Rint Stats] = regress(precip,X);

% T, F tests
% test whether precip has a mean of 34 at the 5% confidence level
[h, p, sig] = ttest(precip,34,0.05);
% test whether precip=month at the 5% confidence level
[h, p, ci] = ttest2(precip,month,0.05);


% Unconstrained Optimization
x0 = [-1,1];
options = optimset('LargeScale','off');
[x,fval,exitflag,output] = fminunc('objfun', x0, options);


% Constrained Optimization
x0 = [-1, 1];
options = optimset('LargeScale', 'off');
[xc,fvalc] = fmincon('objfun', x0, [],[],[],[],[],[], 'confun',options);

% MLE
% The following returns an MLE and a 95% confidence interval for the success probability of a binomial distribution with 20 trials:
% Simulated data, p = 0.75
data = binornd(20,0.75,100,1); 
[phat,pci] = mle(data,'distribution','binomial','alpha',.05,'ntrials',20);
% How to do inference? Need variance/covariance matrix - type help mle

% Other types of fits: e.g. NLS see help stats


% Programming - can only learn by doing
    % try the following
        % for/while loops
        % if/else constructs
        % value function iteration
    % advice: write program by hand before jumping on the computer 
    
% Try to compare stata with matlab - good check