% Script 1 for Week 11

% The purpose of this section is to present a general problem and provide a tool for framing a solution and testing the quality of that solution. 

% The basic idea is as follows. You are collecting information, updated in real-time, about a physical thing. For instance we might have radar information about a plane, or position information about a robot leg. 

% This information will not be perfect. There will be random or systematic errors in the sensory apparatus, limits to your ability to ``read the dial'' and so on.

% Our specific goal here will be to take a finite sequence of data values, assumed to be taken at equally-spaced time intervals, and predict the next value. We will then want to measure how close we got to the ACTUAL next (noisy) value.

% All of this will, of course, depend on the standard deviation and nature of the noise and the specific predictor used.

% We present a variety of predictors below, and all of them take a specific length of data-vector and predict the next value based on the evidence of that data vector. Some make derivative estimates and estimate the subsequent value using that. Others fit the data points to a polynomial and predict using the best-fit polynomial. 

% PredictorPlot.m is our main tool, and can use any of the predictor functions listed below THAT. PredictorPlot.m measures the success of the prediction method on a specific list of noisy data which you provide. This list must be at least as long as the data required by the predictor. PredictorPlot.m  could also use ANY predictor that you might produce which takes a specific length of data-vector and predicts the next value.

% We are trying to get MATLAB to predict the next point in a fairly long noise-added path of a function using a sliding window of ``measured'' data. I test below against values of a known function with normally distributed noise added.  

% PredictorPlot.m does assume you created the data vector using fun, but could be easily modified to use measured data only. In that case the residuals (the errors) would be measured relative to the supplied data vector alone rather than the underlying exact function values. You would create that modified version after you had tested your favorite predictor against all the noisy function-shapes you are likely to encounter and were satisfied by its performance characteristics.

% All of the listed predictors work very very well when there is no noise. And all behave pretty poorly when using noisy data.

% Since PredictorPlot.m is a general purpose tool, if you can come up with a better predictor, PredictorPlot.m will evaluate how well it works on your data.

% Kalman Filtering and PID controllers are subjects many Engineering students learn about and which are used in practice in robotics and signal processing in general. Packages are available with MATLAB to implement or assist with these. IF ONLY I KNEW SOMETHING ABOUT THEM ... I might be able to do better with my predictors! Can you?


clear
t=linspace(1, pi, 10);
f=@(t) 3*sin(t/3)+t.^2;
y=f(t);

z=y+0.3*randn(1,10);

figure(7) 
plot(t,y,t,z,'*');
legend('exact', 'noisy')

PredictorPlot(@(t) f(t), 1, pi, z, @(h) PredLinRegress4(h), 4)

PredictorPlot(@(t) f(t), 1, pi, z, @(h) PredQuadRegress4(h), 4)

PredictorPlot(@(t) f(t), 1, pi, z, @(h) PredCubeRegress5(h), 5)

PredictorPlot(@(t) f(t), 1, pi, z, @(h) Predict4(h), 4)

PredictorPlot(@(t) f(t), 1, pi, z, @(h) Predict3(h), 3)


	
%*********************************%**********************************
function [w,R,Last]=PredictorPlot(fun,a,b,z,predictor,n, stdev)

%  This function takes a list of measured values z at UNIFORMLY spaced times and makes a prediction w of future function values. 

% fun is the function being sampled at regular time intervals on the interval [a,b] at domain values t=a:delta:b where delta=(b-a)/(length(z)-1).

% PredictorPlot.m does assume you created the data vector using fun, but could be easily modified to use measured data only. In that case the residuals (the errors) would be measured relative to the supplied data vector alone rather than the underlying exact function values. You would create that modified version after you had tested your favorite predictor against all the noisy function-shapes you are likely to encounter and were satisfied by its performance characteristics.

% z is a vector of data: the function values fun(t) plus noise (added by the user).

% stdev is optional, used to help make a nicer residuals graph, and is the standard deviation of the noise you added.

% n is the number of consecutive data points required by the specified predictor function.

% The purpose of this function is to test how well a specified predictor can see past the noise to produce future function values of function fun. 

% w is the vector of predicted values and length(w)=length(z)-n so the first entry in w corresponds to time t(n+1).

% R is the residual vector R=w-fun(t(n+1:length(z))).

% A perfect prediction would produce the added noise alone as residual. 

% Last is the prediction for the unmeasured f(b+delta).

% invoke as PredictorPlot(@(t) sin(t/20), 1, pi, z, @(h) Predict4(h), 4)

if nargin<6,error('at least 6 input arguments required'),end

if nargin<7,stdev=0.03;end

if length(z)<n+1 error('The number of data points must exceed the number required by the predictor.'); 
end

% time interval
delta=(b-a)/(length(z)-1);
% data times
t=a:delta:b;
% true function values
y=fun(t);
% length of data vector
L=length(z);
% length of prediction vector
L2=L-n;
% times corresponding to the prediction vector values
t2=t((n+1):L);

for i=(n+1):L;
w(i-n)=predictor(z((i-n):(i-1)));
end

Last=predictor(z((L-n+1):L));

R=w-fun(t2);

figure(58)
plot(t2,R, t2, stdev*ones(length(t2)), t2, -stdev*ones(length(t2)));
legend('residual','+ sigma','- sigma')


figure(59)
clf

plot(t,z,'*r','markersize',12)
hold on
plot(t,y)
plot(t2, w)
legend('noisy data','function','prediction')
hold off

end

%*********************************%**********************************
function pred=Predict2(y)

% This function tries to predict the next value of changing numbers 
% assumed to be a function of a parameter such as time. 
% The time increment between consecutive data points are assumed to be constant.
% This predictor requires two data points and predicts a third.


% Uses a backward difference approximation to the derivative.

if length(y)~=2 error('To use this predictor the data vector must be of length  2') ; end



pred=y(2)+(y(2)-y(1));

end

%*********************************%**********************************
function pred=Predict3(y)

% This function tries to predict the next value of changing numbers 
% assumed to be a function of a parameter such as time. 
% The time increment between consecutive data points are assumed to be constant.
% This predictor requires three data points and predicts a fourth.

% Uses a backward difference approximation to the derivative at the 
% last data point to project the next value.

if length(y)~=3 error('To use this predictor the data vector must be of length 3'); 
end


pred=y(3)+(y(1)-4*y(2)+3*y(3))/2;

end
%*********************************%**********************************
function pred=Predict4(y)

% This function tries to predict the next value of changing numbers 
% assumed to be a function of a parameter such as time. 
% The time increment between consecutive data points are assumed to be constant.
% This predictor requires four data points and predicts a fifth.

% Uses a backward difference approximation to the derivative.

if length(y)~=4 error('To use this predictor the data vector must be of length 4'); 
end

pred=y(4)+(-2*y(1)+9*y(2)-18*y(3)+11*y(4))/6;

end
%*********************************%**********************************
function pred=Predict5(y)

% This function tries to predict the next value of changing numbers 
% assumed to be a function of a parameter such as time. 
% The time increment between consecutive data points are assumed to be constant.
% This predictor requires five data points and predicts a sixth.

% Uses a backward difference approximation to the derivative.

if length(y)~=5 error('To use this predictor the data vector must be of length 5'); 
end

pred=y(5)+(3*y(1)-16*y(2)+36*y(3)-48*y(4)+25*y(5))/12;

end
%*********************************%**********************************
function pred=PredLinRegress4(y)

% This function tries to predict the next value of changing numbers 
% assumed to be a function of a parameter such as time. 
% The time increment between consecutive data points are assumed to be constant.
% This predictor requires 4 data points and predicts a fifth.
% It uses linear least-squares regression.

if length(y)~=4 error('To use this predictor the data vector must be of length  4') ; end


A=polyfit(1:4,y,1);
pred=polyval(A,5);

end
%*********************************%**********************************
function pred=PredCubeRegress5(y)

% This function tries to predict the next value of changing numbers 
% assumed to be a function of a parameter such as time. 
% The time increment between consecutive data points are assumed to be constant.
% This predictor requires 4 data points and predicts a fifth.
% It uses cubic least-squares regression.

if length(y)~=5 error('To use this predictor the data vector must be of length  4') ; end


A=polyfit(1:5,y,3);
pred=polyval(A,6);

end
%*********************************%**********************************
function pred=PredQuadRegress4(y)

% This function tries to predict the next value of changing numbers 
% assumed to be a function of a parameter such as time. 
% The time increment between consecutive data points are assumed to be constant.
% This predictor requires 4 data points and predicts a fifth.
% It uses quadratic least-squares regression.

if length(y)~=4 error('To use this predictor the data vector must be of length  4') ; end


A=polyfit(1:4,y,2);
pred=polyval(A,5);

end
%*********************************%**********************************