function simple_linear_regression
%built by keith regner 10/2016 with R2014b (8.4 OSX)
clear all
close all

%% description
% This demo shows the affect of noise and outliers in linear regression. 
% The line y = x + 1 is shown in the top plot (dashed black line). 
% Experimental data are generated from this equation by adding white 
% Gaussian noise and are shown on the same plot (blue circles). These data 
% are fit with a linear model (solid green line). The residuals of that 
% fitting are shown in the lower plot. Change the amount of ?noise? in the 
% data using the slider bar and observe the change in the equation of best 
% fit, the r2 value, and the residuals plot. To add an outlier, press the 
% space bar and select any location on the top plot to add the point. The 
% point will replace the middle point of the original data set. The fit to 
% the original data and the data with the outlier will both be shown in the 
% upper plot. The residuals of the original and new fits are shown in the 
% lower plot. Observe the change in change in the equation of best fit, the 
% r2 value, and the residuals plot.

%% create data
b.x_data = 1:.25:10;
b.y_data = 2:.25:11;
%add noise to data
b.snr = log10(100);
b.data = [b.x_data' awgn(b.y_data,(10^b.snr))'];
%x axis range over which to fit data         
b.x = 0:.01:12;
slide_min = log10(1);
slide_max = log10(100);

%% fit data
[a] = simp_lin_reg(b.data,b.x);

%% figure and axis formatting
b.fig = figure('Units','inches','Position',[5 .25 8 9],...
    'Name','simple_linear_regression','NumberTitle','off');
b.ax1 = axes('Units','inches','Position',[1 5 6 3.5],'FontSize',12,...
    'LineWidth',1.5,'XLim', [b.x(1) b.x(length(b.x))],'YLim',...
    [0 4+max(b.data(:,2))],'Box','on');
hold on
xlabel('$x$','Interpreter','latex','FontSize',14);
ylabel('$y$','Interpreter','latex','FontSize',14);

b.ax2 = axes('Units','inches','Position',[1 1.3 6 3],'FontSize',12,...
    'LineWidth',1.5,'XLim', [b.x(1) b.x(length(b.x))],...
    'YLim',[-4 4],'Box','on');
hold on
xlabel('$\hat{y}_{i}$','Interpreter','latex','FontSize',14);
ylabel('$y_{i}$ - $\hat{y}_{i}$','Interpreter','latex','FontSize',14);

%% plot solution
%data and fit
b.h_actual = plot(b.ax1,b.x,a.y_hat,'Color',[0 0 0],...
    'LineWidth',.5,'LineStyle','--');
b.h_data = plot(b.ax1,b.data(:,1),b.data(:,2),'o','MarkerFaceColor',...
    [0 0.44 0.74],'MarkerEdgeColor','k','MarkerSize',6);
b.h_fit = plot(b.ax1,b.x,a.y_hat,'Color',[0 1 0],...
    'LineWidth',1.5);
b.h_fit_temp = plot(b.ax1,b.x,a.y_hat,'Color',[0.5 0 0.5],...
    'LineWidth',1.5,'Visible','off');
b.point_old = plot(b.ax1,4,4,'o','MarkerFaceColor',...
    [1 1 0],'MarkerEdgeColor','k','MarkerSize',6,'Visible','off');
b.point_new = plot(b.ax1,4,4,'o','MarkerFaceColor',...
    [1 0 0],'MarkerEdgeColor','k','MarkerSize',6,'Visible','off');
b.point_leg = plot(b.ax1,4,4,'o','MarkerFaceColor',...
    [1 1 1],'MarkerEdgeColor','k','MarkerSize',6,'Visible','off');
b.leg1 = legend(b.ax1,[b.point_leg],{'data'},...
    'Location','NorthWest','Box','off','FontSize',12);

%residuals
b.h_res_data = plot(b.ax2,b.data(:,2),a.res,'o','MarkerFaceColor',...
    [0 0.44 0.74],'MarkerEdgeColor','k','MarkerSize',6);
b.h_res_outlier = plot(b.ax2,b.data(:,2),a.res,'o','MarkerFaceColor',...
    [1 0 0],'MarkerEdgeColor','k','MarkerSize',6,'Visible','off');
b.h_res_ref = plot(b.ax2,b.x,b.x*0,'Color',[0 0 0],...
    'LineWidth',1.5);
b.leg2 = legend(b.ax2,[b.h_res_data b.h_res_outlier],...
    {'without outlier','with outlier'},...
    'Location','NorthWest','Box','off','FontSize',12);
set(b.leg2,'Visible','off');

%% more legend
b.line_act = annotation('line','Units','inches','x',[1.2 1.63],...
    'y',[8 8],'LineWidth',.5,'LineStyle','--','Color',[0 0 0]);
b.txtboxa = annotation('textbox','Units','inches','Position',...
    [1.61 7.8 4 .3],'String','actual','FontSize',12,'LineStyle',...
    'none');

b.line_fit = annotation('line','Units','inches','x',[1.2 1.63],...
    'y',[7.7 7.7],'LineWidth',1.5,'Color',[0 1 0]);
b.txtboxf1 = annotation('textbox','Units','inches','Position',...
    [1.61 7.5 4 .3],'String','fit, ','FontSize',12,'LineStyle',...
    'none');
b.txtbox1 = annotation('textbox','Units','inches','Position',...
    [1.83 7.55 4 .3],'String',...
    ['$\hat{y}$ = ' num2str(a.B_hat1) '\textit{x} + ' ...
    num2str(a.B_hat0) ', $r^{2}$ = ' num2str(a.r^2)],...
    'FontSize',14,'LineStyle','none','Interpreter','latex');

b.line_outlier = annotation('line','Units','inches','x',[1.2 1.63],...
    'y',[7.4 7.4],'LineWidth',1.5,'Color',[0.5 0 0.5],'Visible','off');
b.txtboxf2 = annotation('textbox','Units','inches','Position',...
    [1.61 7.2 4 .3],'String','fit, ','FontSize',12,'LineStyle',...
    'none','Visible','off');
b.txtbox2 = annotation('textbox','Units','inches','Position',...
    [1.83 7.25 4 .3],'String',...
    ['$\hat{y}$ = ' num2str(a.B_hat1) '\textit{x} + ' ...
    num2str(a.B_hat0) ', $r^{2}$ = ' num2str(a.r^2)],...
    'FontSize',14,'LineStyle','none','Interpreter','latex','Visible','off');

%% add user control
b.txtbox = annotation('textbox','Units','inches','Position',...
    [3.7 0.1 4 .3],'String',...
    'noise','FontSize',12,'LineStyle','none');
b.sl = uicontrol('Parent',b.fig,'Style','slide','Units','inches',...
    'Position',[1.5 .5 5 .2],'Min',slide_min,'Max',slide_max,'Value',...
    b.snr,'SliderStep',[1/500 1/50],'Tag','slider1',...
    'BackgroundColor',b.fig.Color,'UserData',...
    struct('DATA',b.data,'X',b.x),'Callback',{@slider_call,b});
b.label1 = uicontrol('Parent',b.fig,'Units','inches','Style',...
    'text','Position',[.8 .4 .6 .3],'String','a lot',...
    'FontSize',12,'BackgroundColor',b.fig.Color);
b.label2 = uicontrol('Parent',b.fig,'Units','inches','Style',...
    'text','Position',[6.5 .4 .8 .3],'String','a little',...
    'FontSize',12,'BackgroundColor',b.fig.Color);

set(b.fig,'WindowKeyPressFcn',{@pb_kpf,b,b});

set(findall(b.fig,'-property','FontUnits'),'FontUnits','normalized')
set(findall(b.fig,'-property','Units'),'Units','normalized')


end

%% callback functions
%window key press
function [] = pb_kpf(varargin)
[b,c] = varargin{[3,4]};
old_data = c.data;
data = c.data;
[xx,yy] = ginput(1);
oldx = data(round(length(data(:,1))/2),1);
oldy = data(round(length(data(:,2))/2),2);
data(round(length(data(:,1))/2),1) = xx;
data(round(length(data(:,2))/2),2) = yy;
new_data = [data(:,1) data(:,2)];
[c] = simp_lin_reg(new_data,b.x);
b.old_data = old_data;

set(b.point_old,'xdata',oldx,'ydata',oldy,'Visible','on');
set(b.point_new,'xdata',xx,'ydata',yy,'Visible','on');
set(b.h_data,'xdata',new_data(:,1),'ydata',new_data(:,2))
set(b.h_fit_temp,'ydata',c.y_hat,'Visible','on')
set(b.h_res_outlier,'xdata',new_data(:,2),'ydata',c.res,...
    'Visible','on')
set(b.h_res_data,'MarkerFaceColor',[1 1 0])
set(b.leg2,'Visible','on');
set(b.line_outlier,'Visible','on')
set(b.txtboxf2,'Visible','on')
set(b.txtbox2,'Visible','on','String',...
    ['$\hat{y}$ = ' num2str(c.B_hat1) '\textit{x} + ' ...
    num2str(c.B_hat0) ', $r^{2}$ = ' num2str(c.r^2)])
end

function slider_call(varargin)
%callback for the slider
[H,b] = varargin{[1,3]};
slider_value = (10^(H.Value));
set(b.point_old,'Visible','off')
set(b.point_new,'Visible','off')
set(b.h_fit_temp,'Visible','off')
set(b.line_outlier,'Visible','off')
set(b.txtboxf2,'Visible','off')
set(b.txtbox2,'Visible','off')
set(b.leg2,'Visible','off');
set(b.h_res_outlier,'Visible','off')
set(b.h_res_data,'MarkerFaceColor',[0 0.44 0.74])
data = [b.x_data' awgn(b.y_data,slider_value)'];
[c] = simp_lin_reg(data,b.x);
c.data = data;

set(b.h_data,'xdata',data(:,1),'ydata',data(:,2))
set(b.h_fit,'ydata',c.y_hat)
set(b.h_res_data,'xdata',data(:,2),'ydata',c.res)
set(b.txtbox1,'String',...
    ['$\hat{y}$ = ' num2str(c.B_hat1) '\textit{x} + ' ...
    num2str(c.B_hat0) ', $r^{2}$ = ' num2str(c.r^2)])
set(b.fig,'WindowKeyPressFcn',{@pb_kpf,b,c});

end

%% simple linear regression function
function [a] = simp_lin_reg(data,x)
a.x_bar = mean(data(:,1));
a.y_bar = mean(data(:,2));
a.B_hat1 = sum((data(:,1)-a.x_bar).*(data(:,2)-a.y_bar))/...
    sum((data(:,1)-a.x_bar).^2);
a.B_hat0 = a.y_bar - a.B_hat1*a.x_bar;
a.y_hat = a.B_hat0 + a.B_hat1*x;
a.r = sum((data(:,1)-a.x_bar).*(data(:,2)-a.y_bar))/...
    (sqrt(sum(((data(:,1)-a.x_bar).^2)))*...
    sqrt(sum(((data(:,2)-a.y_bar).^2))));
y_hat_res = a.B_hat0 + a.B_hat1*data(:,1);
a.res = data(:,2) - y_hat_res;
end

