function [ h, lms ] = plot_regression( x, y, xname, yname, clr, mrk_size, xlims, filled )
% plot_regression - Plots a regression line with confidence interval for 
% vectors x and y

if ~iscell(x)
   x = {x};
   y = {y};
end

if nargin < 3
    xname = 'x';
end

if nargin < 4
    yname = 'y';
end

if nargin < 5
   clr = lines(size(x,2));
end

if nargin < 6
   mrk_size = 70; 
end
    
if nargin < 7
    minx = Inf;
    maxx = -Inf;

    for i = 1 : size(x,2)
        x_i = x{i};
        minx = min(minx, min(x_i));
        maxx = max(maxx, max(x_i));
    end

    minx = minx - minx*.05;
    maxx = maxx + maxx*.05;
else
    minx = xlims(1);
    maxx = xlims(2);
end

if nargin < 8
   filled = true; 
end

dx = (maxx - minx) / 100;
xx = (minx:dx:maxx)';

lms = cell(0,size(x,2));

for i = 1 : size(x,2)

    x_i = x{i};
    y_i = y{i};
    lclr = lighten_colour(clr(i,:), 0.5);

    data = array2table([x_i,y_i],'VariableNames',[{xname},{yname}]);
    model = sprintf('%s~%s', yname, xname);

    lm = fitlm(data, model);
    lms(i) = {lm};
    
    [ypred, ci] = predict(lm, xx);
    hh = fill([xx;flip(xx)], [ci(:,1);flip(ci(:,2))], lclr);
    hh.LineStyle = 'none';
    hh.FaceAlpha = 0.5;
    hold on;
    
end

for i = 1 : size(x,2)
    if filled
        h = scatter(x{i}, y{i}, mrk_size, clr(i,:), 'filled');
    else
        h = scatter(x{i}, y{i}, mrk_size, clr(i,:));
    end
end

for i = 1 : size(x,2)
    ypred = predict(lms{i}, xx);
    hh = line([minx maxx],[ypred(1) ypred(end)]);
    hh.Color = clr(i,:);
    hh.LineWidth = 2;
end

xlim([minx maxx]);

if nargin > 2
    hh = xlabel(xname);
    hh.FontSize = 14;
end

if nargin > 3
    hh = ylabel(yname);
    hh.FontSize = 14;
end

end

