%Lab 9: Resampling methods: Bootstrapping, Monte Carlo, Overfitting

%% 0 Init
clear all
close all
clc

%Initialize the random number generator with a fine-grained system clock
%seed
s = RandStream('mt19937ar','Seed',sum(round(clock*1000))); 
RandStream.setGlobalStream(s); 

%% 1 Classic Efron - does Aspirin prevent heart attacks? 
nA = 11037; %Number of aspirin takers
nHA = 104; %Number of aspirin takers with a heart attack
nP = 11034; %Number of placebo takers
nHP = 189; %Of those, number of heart attacks

nSamples = 1e3; %Number of times to resample - "large-ish" 
ratio = zeros(nSamples,1); %Preallocate memory - makes it faster
widthCI = 95; %How wide should the confidence interval be? 

%Now: a) Make or "represent" the data vectors
A1 = ones(nHA,1); %Aspirin takers with heart attacks
A2 = zeros(nA-nHA,1); %Aspirin takers without HA
P1 = ones(nHP,1); %Placebo takers
P2 = zeros(nP-nHP,1); %Same thing

A = [A1;A2]; %Stack them
P = [P1;P2]; %Concatenate

%b) Calculate the distribution of sample means by sampling with replacement
for ii = 1:nSamples %We will do this without a loop soon
    indicesA = randi(nA,[nA 1]); %Draw random integers from 1 to n, with the same sample size as the original sample
    indicesP = randi(nP,[nP 1]); %Draw random integers from 1 to n, with the same sample size as the original sample
    ratio(ii) = sum(A(indicesA))/sum(P(indicesP));
end

%Alternative way - doing this without actually making fake data vectors in
%a loop - which takes a lot of time. This is much faster
A = rand(nA,nSamples); %Create random numbers from 0 to 1, nA timesxnSamples times
P = rand(nP,nSamples); %Same

Atemp = A < nHA/nA; %If it is smaller, render it 1
Ptemp = P < nHP/nP; %Same
ratio = sum(Atemp)./sum(Ptemp); 

%c) Now that we bootstrapped the distribution of ratios, we can calculate the CI
sortedRatio = sort(ratio); %Sort the vector to get a monotonously rising discrete distribution
lowerBoundIndex = round(nSamples/100*(100-widthCI)/2); %Index of the lower bound
upperBoundIndex = round(nSamples/100*(100-((100-widthCI)/2))); %Index of the lower bound

lowerBound = sortedRatio(lowerBoundIndex); %Find actual lower bound value
upperBound = sortedRatio(upperBoundIndex); %Same for upper

%d) Plot histogram
figure
nBins = 50;
histogram(ratio,nBins) %However many bins
xlabel('Ratio')
line([mean(ratio) mean(ratio)], [min(ylim) max(ylim)], 'color','r','linewidth',2)
line([lowerBound lowerBound], [min(ylim) max(ylim)], 'color','g','linewidth',2)
line([upperBound upperBound], [min(ylim) max(ylim)], 'color','g','linewidth',2)
title('Bootstrapped distribution of heart attack ratios Aspirin vs. Placebo')

%% 2 Using resampling methods for cross-validation
%a) Get the samples
%b) Fit polynomials
%c) Leave one out and resample

numPoints = 8; %Eero's choice
noiseMagnitude = 2; %Some dial that allows us to tune up and down the noise
x = -5:0.01:5; %1001 equally spaced points between -5 and 5
y = x.^2 + noiseMagnitude .* randn(1,length(x)) ;%Ground truth: Parabola, but with additive random noise

figure
plot(x,y)
title('Ground truth: A noisy parabola')

%a) Find the 8 points on this parabola
samplingIndices = randi(length(x),[numPoints,1]); %Find the indices of the sampling points
hold on;
plot(x(samplingIndices),y(samplingIndices),'.','color','r','markersize',20)
shg

%b) (Over)fitting polynomials and calculating RMSE - which will be
%misleading
RMSE = []; %Our figure of merit
figure
for ii = 1:numPoints
    subplot(2,4,ii)
    p = polyfit(x(samplingIndices),y(samplingIndices),ii); %Use randomly drawn data points and fit it with polynomials of degree ii
    %p will be a vector of weights/coefficients of the polynomial
    yHat = polyval(p,x); %Evaluate these estimated p's at all values of x
    plot(x,yHat)
    hold on
    plot(x(samplingIndices),y(samplingIndices),'.','color','r','markersize',20)
    title(num2str(ii))
    RMSE(ii) = sqrt(mean((y(samplingIndices)-yHat(samplingIndices)).^2)); 
    %RMSE does not represent the sign of the error. If that matters, don't
    %use RMSE
end

%This seems to work very well, by visual inspection. Quantifying it with
%RMSE: 
figure
plot(1:numPoints,RMSE,'linewidth',2)
title('RMSE as a function of degree of polynomial')
xlabel('Degree of polynomial')
ylabel('RMSE')

%As you can see, the RMSE keeps dropping as we add terms to the polynomial.
%This is necessarily the case, but we are fitting to noise. Overfitting.
%So this won't generalize

%Solution: Cross-validation.
%Big idea: Split data into a training set to build your model and a test
%set to test it. How to split it is up for debate. Usually 80/20, 90/10, or
%"leave one out". 

RMSE = zeros(nSamples,numPoints-1); %Reinitialize RMSE

for ii = 1:nSamples %Do this as many times as you can afford, computationally
    %Splitting into test and training set. Many issues: Which to do first,
    %how to split - in sequence or random. That's all fungible. As long as
    %you do it often enough
    testIndex = randi(numPoints,1); %Index of the point we want to predict
    testSet = samplingIndices(testIndex); %Indices of the test-set
    trainingSet = samplingIndices; %Make a copy of the samplingIndices
    trainingSet(testIndex,:) = []; %Get rid of it - We could do this in many other ways
    for jj = 1:numPoints-1
        p = polyfit(x(trainingSet),y(trainingSet),jj); %Fit to training set, not the full set(!)
        yHat = polyval(p,x); %Evaluate anew
        RMSE(ii,jj) = sqrt((y(testSet)-yHat(testSet)).^2); 
    end
end

%LOGIC/INSIGHT: FIT to the training set, but calculate the ERROR with the
%test set
figure
plot(1:numPoints-1,mean(RMSE))
title('REAL RMSE as a function of degree of polynomial')
xlabel('Degree of polynomial')
ylabel('RMSE')

%% 3 Permutation tests to compare sample means as a replacement for significance tests
%a) Enter the data and pick parameters, e.g. number of times you want to do
%this
%This is our rat example from the lecture

K = [117 123 111 101 121]; %Number of dendritic spines in Ketamine treated neurons
C = [98 104 106 92 98]; %Number of dendritic spines per control neuron

numMCS = 1e4; %Number of "Monte Carlo" simulations
nK = length(K); %Number of Ketamine treated neurons
nC = length(C); %Number of control neurons

%b) Come up with a reasonable test statistic
testStatistic = sum(K) - sum(C); 
%Thisis something that requires thought and creativity. Depends on your
%hypothesis. This could be one-tailed or two-tailed or whatever
%The direction of the difference might matter. 
%Also, say your hypothesis pertains to whether Ketamine multiplies the
%number of dendrites, vs. adding some, a *ratio* might be the more
%appropriate test statistic. 

%Our empirical test statistic is 75. 
%We now need to assess how likely that is by chance
%Big idea: We permute the values and form arbitrary groups (pretending we
%didn't know which neuron was in the treatment group vs. control group)
%Doing this a large number of times will allow us to create a "null
%distribution" 
%Logic: If Ketamine had no effect, it should not matter which group the
%neurons came from 

%c) Do the Monte Carlo method
dTS = zeros(numMCS,1); %Initialize and preallocate the vector that will represent the distribution of teststatistics
allData = [K C]; %Concatenate all the data

%Monte Carlo method: Random permutation/reordering/rearrangement of the
%group labels

for ii = 1:numMCS %However many times you want
    indices = randperm(nK+nC); %Reshuffle as many times as there are total datapoints
    %Bootstrap: Sampling with replacement
    %MC: Rearranging or sampling without replacement
    tempDATA = allData(indices); %Create a newly permuted dataset
    %Recarving two groups
    group1 = tempDATA(1:nK); %Capture the first n elements, our virtual experimental group
    tempDATA(:,1:nK) = []; %Delete those from teh permuted dataset
    group2 = tempDATA; %Group 2 is simply what is left
    dTS(ii) = sum(group1)-sum(group2) ; %Capture the simulated test statistic 
end

figure
histogram(dTS,nBins)

%d) Find the *exact* p value
%By assessing the proportion of times the randomly shuffled teststatistic
%exceeded the empirical one

exactP = sum((dTS >= testStatistic))/length(dTS); %Proportion, this is one-tailed
exactP = sum(testStatistic<=abs(dTS))/length(dTS); %Two tailed case

line([testStatistic testStatistic],[min(ylim) max(ylim)],'color','k','linewidth',3)
xlabel('Distribution of simulated test statistics in absence of an effect')
title(['Distribution of test statistics. Exact p = ', num2str(exactP,'%1.3f')])
set(gcf,'color','w')

%% 4 To add: Stability of means, manually and with in-built functions.

