% This livescript will cover the topics of Cross Validation and Bootstrap Respampling.  Both of these methods are ways of using a single dataset (cleverly) to get more information about the fit that you make to that data than might otherwise be available; e.g., among many methods, which would have the lowest test error (recall the bias variance tradeoff), and how confident am I in the actual fit that is being produced?
% Let's start with Cross Validation first.  Its main goal is to use a single dataset in a particular way in order to estimate the test error as a function of model flexibility (for example), to aid in model selection.  Before we get to Cross Validation and how it works, recall from the discussion and livescript about the bias variance tradeoff that increasing the flexibility of our model will generally always decrease training error, but that test error typically exhibits a "U" shaped curve as a function of model flexibility, with a (generally unknown) level of flexibility existing that minimizes this test error.
% Borrowing some code from the bias variance livescript, we can illustrate this behavior:

xa=0; %this is the left hand side of our x domain
xb=9; %this is the right hand side of our x domain
n=50; % the number of data points
sig2=.5;
f=@(x) 4+15*x.*exp(-x); %this is just an example - try some other function if you want!
xvals0=linspace(xa,xb,n)';
ytrue=f(xvals0);
epsilon=normrnd(0,sqrt(sig2),n,1);
yvals0=ytrue+epsilon;

trainmses=zeros(100,n-3); %storage for our trainmse values
testmses=trainmses; %storage for our testmse values
for trial=1:100 %we'll do 100 randomly generated training sets, all using the same f and x values
    epsilon=normrnd(0,sqrt(sig2),n,1);
    yvalst=ytrue+epsilon;
    for order=1:n-3 %for each training set, look at all possible orders of fit
        if order<=3
            spln=spap2(1,order+1,xvals0,yvalst);
        else
            spln=spap2(order-2,4,xvals0,yvalst);
        end
        splntrainerr=fnval(spln,xvals0)-yvalst;
        trainmses(trial,order)=splntrainerr'*splntrainerr/n;
        splntesterr=fnval(spln,xvals0)-yvals0; %here we just use the initial dataset generated in this livescript as the test data
        testmses(trial,order)=splntesterr'*splntesterr/n;
    end
end
xord=1:n-3; %just setting the x values for our plot, which are just the orders used
plot(xord,mean(trainmses),'bo',xord,mean(testmses),'ro',[1,n-3],[sig2,sig2]);
[~,optorder]=min(mean(testmses))

% After running the above code, we find that the optimal order of our spline fit is given by the variable "optorder".  This was found by generating a lot of test datasets from the known function f, then fitting to each and evaluating the test MSE on the test dataset.  We could make our estimate of the optimal order as accurate as we want for this function (on this interval of x values) by using an arbitrarily large number of training and test datasets, since we can always just create more; the only limitation is our patience and computing power!
% But in a real-world problem, we will typically just have one (hopefully large) dataset to use, and we will have no idea what the true function is, so we cannot just generate more datasets.  In theory we could maybe generate more, new datasets by running whatever "experiment" generated the initial dataset again and again, but this might be prohivitively costly or perhaps not even possible (the data isn't from an experiment per se, and cannot really be replicated again even if we had the resources).
%
% So, what are we to do if we want to determine the optimal fitting method for our dataset?  This is where Cross Validation comes in.  Note that there is no reason why we absolutely must use our full datset for fitting our estimate of the function f.  Of course, there are good reasons to do so, namely that more data generally always leads to a better fit, so we naturally want to use as much as possible for training purposes.  But, life is full of tradeoffs, and if we want to estimate test error using a single dataset, we will be forced to reserve some part of the dataset for testing (validation), using what remains for training.
% Cross Validation does precisely this, but in a systematic way.  A specific implementation of Cross Validation will start by choosing an integer k between 2 and n, where n is the total number of datapoints.  The dataset will then be separated (typically randomly) into k non-overlapping subsets.  So, for k=2, the dataset will be cut in half, leaving us with two distinct datasets (but each presumably generated by the same f); for k=n, each individual point becomes its own tiny dataset.  Then, each resulting subset (of which there are k) is used one at a time as a distinct test dataset, where in each of these test cases all other data (the remaining k-1 subsets) are used to train the model.  After all k subsets have been used to test, the average of their test MSE can be used to estimate the true test MSE for the model employed.  One could do this for a variety of models/methods (or levels of flexibility), then pick the one with the lowest average test MSE as the "best".
% Let's try this for the test dataset created by the code above:

k=10; %specify the k value for CV.  Typical choices are 5 or 10.
shuffle=randperm(n); %this will be used to place our data in a random order
xvals=xvals0(shuffle); %randomly permute x vals
yvals=yvals0(shuffle); %randomly permute y vals same as x vals
newn=floor(n/k)*k; %to make things clean, we may drop some datapoints so that n is an integer multiple of k
n=newn;
xvals=xvals(1:n);
yvals=yvals(1:n);
ntest=n/k; %the size of each test set
ntrain=ntest*(k-1); %the size of each training set
xCV=reshape(xvals,k,ntest); %in xCV, each row is one of the k subsets
yCV=reshape(yvals,k,ntest); %ditto for yCV

testmsesCV=zeros(k,ntrain-3); %storage for testmses
for ind=1:k %loop over all subsets as test set
    xtrain=reshape(xCV(setdiff(1:k,ind),:),ntrain,1); %the x training set this time
    ytrain=reshape(yCV(setdiff(1:k,ind),:),ntrain,1); %the y training set this time
    xtest=xCV(ind,:)'; %the x test set this time
    ytest=yCV(ind,:)'; %the y test set this time
    for order=1:ntrain-3 %for each training set, look at all possible orders of fit
        if order<=3
            spln=spap2(1,order+1,xtrain,ytrain);
        else
            spln=spap2(order-2,4,xtrain,ytrain);
        end
        splntesterr=fnval(spln,xtest)-ytest; %test error
        testmsesCV(ind,order)=splntesterr'*splntesterr/ntest;
    end
end
xordCV=1:ntrain-3; %just setting the x values for our plot, which are just the orders used
plot(xord,mean(testmses),'ro',xordCV,mean(testmsesCV),'go',[1,n-3],[sig2,sig2]);
[~,optorderCV]=min(mean(testmses))

% Go ahead and play with the above, choosing different values for k to see what happens.  You could also change the function f above.  For the f that is provided, the optimal order identified by CV is typically quite close to that identified by the simulation method above.
%
% The second topic of this livescript is Bootstrap Resampling.  The main purpose of this method is to answer questions like: how confident am I in the fit that I have estimated from my data?  For some fitting methods, such as ordinary least squares, formulas exist to provide confidence intervals for each of the estimated parameters to help answer this question.  However, not all fitting methods have such readily available formulas.  And, even if these formulas to exist, they may be based on assumptions that don't hold for some particular dataset you are considering.
% Bootstrap Resampling will try to answer this question without any need of specific formulas.  In some ways it is similar to Cross Validation, as Bootstrap Resampling will also consider some number of datasets that are constructed using the initial full dataset, then fit each one of these to the method under consideration.  One big difference is that in Bootstrap Resampling, we actually create an entirely "new" dataset from our full initial dataset, rathern than just taking (random) subsamples like Cross Validation does.  Specifically, in Bootstrap Resampling we will construct B new bootstrap datasets, each of size N (the size of the initial dataset), from the initial dataset.  Each of these B bootstrap datasets is made by randomly selecting, with replacement, N points from the initial dataset.  That is, each point of a bootstrap dataset is made by just randomly selecting (uniformly) one of the points in the initial dataset, and when we generate the next point in the bootstrap dataset we again randomly select (uniformly) one of the points in the initial dataset, even if that point was already selected previously (this is the "with replacement" noted above).  This is how this might look in code:

B=1000; %number of bootstrap samples to construct
sample=randi(n,n,B); %this generates B different random choices of n integers from 1 to n, which is what will indicate which points are in which bootstrap sample;
bootstrapx=xvals0(sample); %our bootstrap x values, with n rows and B columns
bootstrapy=yvals0(sample); %our bootstrap y values, with n rows and B columns

% Before proceeding on to how these bootstrap datasets might be used, it is worth considering for a moment why creating these datasets in the first place might be reasonable.  Suppose you take the view that your original dataset is a set of N random samples from some unknown joint probabilty distribution p_\beta(x,y).  Here, the subscript \beta indicates that this probability distribution likely contains some parameters that have some true values, which we will refer to as \beta^*, but we don't generally know what those true values are.  Some of this relates to what we discussed before about Bayesian regression, so that this interpretation of our dataset is quite reasonable even if we think that the relationship between x and y is a totally deterministic function f, since there will always be the random error \varepsilon at each point.  Under this view, what Bootstrap Resampling is doing is trying to generate new datasets from the unknown distribution p_\beta(x,y) by using the empirical distribution (our original dataset) as a proxy.
% Now, given a particular bootstrap dataset b, we could then attempt to determine the most likely value(s) for the parameter(s) \beta in the unknown probability distribution p_\beta(x,y) by using some fitting method, e.g., OLS.  This will yield some estimate \hat{\beta}_b, which we hope is not so far off from the true values \beta^*.  It is important to recognize that, since the bootstrap (and original!) dataset is recognized to be a set of random variables from some distribution, then the resulting fit of the parameters \hat{\beta}_b must also be a random variable from some other distribution.  That is, we would not expect to get the exact same fit for each bootstrap sample, since the datasets will not all be identical.  Well, it turns out that as N\to\infty (and under some other assumptions we won't detail here), the distribution of our parameter fits is a multinormal distribution:
% \hat{\beta}\sim N[\beta^*,F^{-1}(\beta^*)],
% where F(\beta^*) is the Fisher Information matrix.  Plausibly, then, if we were to construct many bootstrap datasets, fit each one to obtain a value for \hat{\beta}_b, and then look at the distribution of these \hat{\beta}_b values obtained, that distribution's mean may well be close to the true parameter values, and its variance might provide information about how confident we are in that value (and even things like 95% confidence intervals could be obtained, etc).
% In order to explore this, we will perform OLS fits of our bootstrap datasets to the true functional form for this problem to try to estimate the parameters of the true function f (which should be 4 and 15 for the default function).  After obtaining all of these OLS coefficients across our bootstrap samples, we can visualize by making histograms of them, and looking at means, etc.

nformbootstrapx=bootstrapx.*exp(-bootstrapx); %this is the true functional form used in the example here.  There is then a truly linear relationship between bootstrapy and nformbootstrapx.
betas=zeros(2,B); %storage for the various coefficient estimates for out bootstrap samples
for sample=1:B
    xforfit=nformbootstrapx(:,sample);
    xforfit=[ones(n,1),xforfit]; %add our usual column of 1s for the intercept
    yforfit=bootstrapy(:,sample);
    betas(:,sample)=xforfit\yforfit; %Matlab command for least squares solving
end

lower=floor(.05*B)+1; %find index for lowest 5% of values, to use finding 90% CI
upper=floor(.95*B); %find index for highest 5% of values, to use finding 90% CI
histogram(betas(1,:));
xlabel("beta_0");
ylabel("frequency");
fprintf('Mean beta_0 value is %g',mean(betas(1,:)))
ord=sort(betas(1,:));
fprintf('The beta_0 90 percent CI range is %g to %g and standard error is %g',ord(lower),ord(upper),std(ord))
histogram(betas(2,:));
xlabel("beta_1");
ylabel("frequency");
fprintf('Mean beta_1 value is %g',mean(betas(2,:)))
ord=sort(betas(2,:));
fprintf('The beta_1 90 percent CI range is %g to %g and standard error is %g',ord(lower),ord(upper),std(ord))

%
% Just as a point of comparison, we can also run a Matlab built-in linear regression that will report the coefficients and their standard errors, in comparison with those found above:

fitlm(xvals0.*exp(-xvals0),yvals0)