% This script will illustrate the fundamental bias-variance tradeoff in data science through an example.
% Suppose that our variables x and y truly related via some relationship y=f(x).  We obtain a dataset containing n datapoints \left(x_i ,y_i \right), where there is some inherent uncertainty in the measurements such that y_i=f(x_i)+\epsilon_i, where \epsilon_i is a random variable whose distribution has zero mean and a variance of \sigma^2 (we'll just assume a normal distribution for this).
% Below, we set up some of these variables (and a few more) in code for our example:

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!
xvals=linspace(xa,xb,n)';
ytrue=f(xvals);
epsilon=normrnd(0,sqrt(sig2),n,1);
yvals=ytrue+epsilon;

% Now, suppose we were to use some data science method to fit this dataset.  We of course don't know what f is, but we choose some method that we think can generally capture it "well".  This could be a simple linear regression, if we think that's appropriate, or could be some higher order polynomial if we think the true function is very nonlinear.  For the sake of this example, we will choose splines for this fitting task, because it is easily adaptable to range from linear to essentially as nonlinear as we want.  Here is some code that lets us choose the "order" of our spline fit (order 1 is linear, and higher order essentially look like higher order polynomials) then plots what we obtain and computes the Mean Square Error (MSE) of our fit to the data:

order=12;
order=min(order,n-3); %technically there is a maximum order that the method allows, and at that maximum order the fit will pass perfectly through all datapoints
if order<=3
    spln=spap2(1,order+1,xvals,yvals);
else
    spln=spap2(order-2,4,xvals,yvals);
end
%the variable spln is the fit we obtain

xuni=linspace(xa,xb,1000)'; %just for plotting purposes to get a smooth curve
yuni=f(xuni);

figure(1);
plot(xvals,yvals,'bo',xuni,yuni);
hold on;
fnplt(spln);
hold off

splntrainerr=fnval(spln,xvals)-yvals; % the error between our fit and the observed y values
trainmse=splntrainerr'*splntrainerr/n % the mean square error, with lower meaning a better fit

% In this plot, the blue circles are the data points, the red curve is the true function f (which again, we don't know in a real data science problem), and the yellow curve is our fit to the data using our chosen method.  Play around with the "order" above to see how this changes the perceived fit and the MSE.
%
% After playing around with the example above, you should find that in general the MSE gets smaller as you increase the order of your fit.  This makes sense: the higher order makes the fit pass closer and closer to the datapoints, leading to low errors. Clearly, a higher order fit is better able to adapt to a wide variety of functions f.  However, generally speaking, no one is actually too concerned with the MSE between the fit and the training data used to make that fit.  This is precisely because this error can be made as small as you want (including 0) by simply choosing a high enough order method.
% What people really care about is how far away your fit is from the true function f.  Unfortunately, since we never know the true function f, this is not really possible to directly measure.  But, we can get a proxy measure for this by considering what would happen if we compared our fit from above to a second, independent testing dataset that has the same underlying function f.  If our fit from above is actually quite close to the true function, then we would expect that another dataset drawn from that same function would have a relatively small MSE, at least roughly like the MSE between our fit and the training dataset.  Let's try this below:

order=8;
order=min(order,n-3); %technically there is a maximum order that the method allows, and at that maximum order the fit will pass perfectly through all datapoints
if order<=3
    spln=spap2(1,order+1,xvals,yvals);
else
    spln=spap2(order-2,4,xvals,yvals);
end
%the variable spln is the fit we obtain

splntrainerr=fnval(spln,xvals)-yvals; % the error between our fit and the observed y values
trainmse=splntrainerr'*splntrainerr/n % the mean square error, with lower meaning a better fit

%we'll just use the same x values for the test dataset
epsilon=normrnd(0,sqrt(sig2),n,1); % a new set of errors for our test set
ytest=ytrue+epsilon; %our test y values
splntesterr=fnval(spln,xvals)-ytest; % the error between our fit and the test y values, which are not used to construct the fit but are generated by the same f
testmse=splntesterr'*splntesterr/n % the mean square error, with lower meaning a better fit

% Again, play around with the order above.  In any case, you should generally find that the test MSE is greater than the train MSE, which of course makes sense.  However, what is really interesting is the way that the train MSE behaves as we change the order.  We will study this more systematically below.
%
% Let's now study a little more systematically how the trainmse and testmse vary with order of fit.  We'll do so by looping over several possible training datasets and several possible orders, comparing all of these fits to a single test set, then look at averages of the various MSE obtained as a function of order.

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,xvals,yvalst);
        else
            spln=spap2(order-2,4,xvals,yvalst);
        end
        splntrainerr=fnval(spln,xvals)-yvalst;
        trainmses(trial,order)=splntrainerr'*splntrainerr/n;
        splntesterr=fnval(spln,xvals)-yvals; %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
figure(2);
plot(xord,mean(trainmses),'bo',xord,mean(testmses),'ro',[1,n-3],[sig2,sig2]);

% In this plot, we see three "curves".  The blue circles are the training MSE, which just decrease with increasing order, as discussed above.  The red curve is the test MSE, which (for this example) appears to decrease with order at first, but then increases.  The yellow horizontal line is showing the location of \sigma^2: in theory, the expected test MSE could never be less than this, as the data inhenrently has at least this much error in it.
% This plot is a visual representation of the bias-variance tradeoff in data science.  Here, the term bias refers to how well a particular method is able to fit a variety of functions.  A more flexible method (synonymous with higher order here) can fit a complicated nonlinear function as well as a linear one, but a less flexible method (like linear regression) can only fit less complicated functions well.  So, it would seem that choosing a very flexible method would be good, as we can then possibly fit our unknown function even if it is quite nonlinear.  This is best understood by the fact that the training MSE gets smaller and smaller as order increases here.
% However, as illustrated by the test MSE in our plot, this does not tell the whole story.  This is because these highly flexible methods are prone to displaying a lot of variance (the other half of the bias-variance tradeoff).  What variance means here is how different a given fit will be from one dataset to the next, even if the datasets all come from the same underlying f.  As seen above, when we use a very flexible method (high order), the fit passes through all the training datapoints perfectly.  But, note that the training points are not all actually on the true curve f, due to the inherent error \epsilon!  So, trying to fit training data too closely will actually make the fit vary substantially from one training set to the next, leading to high variance and consequently a high(er) test MSE than a lower flexibility method.
% We can then understand the behavior of the red dots in our plot.  Our function f in this example is not very linear, so a low flexibility method will have high bias (not close to the true f), but small variance.  The resulting test MSE is somewhat high then, taking both into account.  On the other hand, a very flexible method will have low bias (can fit our nonlinear function quite well), but high variance (each training set will give a notably different fit), again leading to relatively high test MSE.  Somewhere in between is a level of flexibility that best balances the tradeoff between bias and variance, yielding the smallest test MSE of all, which cannot be any smaller than \sigma^2 in any case.