-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathRegression
110 lines (99 loc) · 2.67 KB
/
Regression
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
clear all;
s=6; % set random number seed.
rand('seed',s); randn('seed',s);
% Make noisy data with slope m and intercept c.
% Define 10 salary values for horizontal axis.
s = 1:1:11; s=s';
% Set true values of slope m and intercept c.
m = 0.5;
c = 3;
% Set standard deviation of each measured height.
sds = 2*[1:11]'/10; %rand(size(s))*2;
sds = sds .* [1:11]'/10;
% Use sds to make vector of noise.
eta = randn(size(s)).*sds;
% Pretend some data points are extreme.
eta(8)=-1; eta(9)=-3; eta(11)=-3;
% Find observed values of x (with noise added).
x = m*s + c + eta;
% Weighted Least Squares regression. We could find the solution using the smallest value in the Farray below also.
% Find weightings w (discount) for each data point.
vars = sds.^2;
w=1./vars;
%w=ones(size(w)); % un-comment this line to get solution based on uniform noise terms.
ss=[ones(size(s)) s]; % prepend column of ones so that solution includes intercept term.
[params,stdx,mse,S] =lscov(ss,x,w);
mest2 = params(2);
cest2 = params(1);
%%%%%%%%
xest2 = mest2.*s + cest2;
c0 = cest2;
m0 = mest2;
% Plot fitted line xest (=xhat in text) and data points.
figure(1); clf;
plot(s,x,'k*',s,xest2, 'k', 'LineWidth',2, 'MarkerSize',10);
set(gca,'Linewidth',2); set(gca,'FontSize',20);
xlabel('Salary, {\it s} (groats)');
ylabel('Height, {\it x} (feet)');
set(gca,'XLim',[0 12],'FontName','Ariel');
set(gca,'YLim',[0 9],'FontName','Ariel');
% Plot sd from each data point
for i=1:11
x1 = x(i)-sds(i); x2=x(i)+sds(i);
s1 = i; s2 = i;
ss = [s1 s2]; xx = [x1 x2];
hold on;
plot(ss,xx,'k','LineWidth',2); hold off;
end
% Code for 2d plot
m = mest2;
mmin = mest2-1; % 0.5
mmax = mest2+1;
c = cest2;
cmin = cest2-1; %0.4
cmax = cest2+1;
minc = (mmax-mmin)/100;
cinc = (cmax-cmin)/100;
Fs = [];
ms = mmin:minc:mmax;
nm = len(ms);
cs = cmin:cinc:cmax;
nc = len(cs);
Farray = zeros(nm,nc);
for m1 = 1:nm
for c1 = 1:nc
mval=ms(m1);
cval=cs(c1);
y1 = mval*s + cval;
F1 = ((x-y1)./sds).^2;
Farray(m1,c1)=sum(F1);
Fs = [Fs F1];
end
end
figure(2);
Z1 = Farray';
zmin=min(Z1(:));
zmax=max(Z1(:));
range=zmax-zmin;
v=zmin:range/10:zmax;
% Adjust spacing of contour lines.
v=0:0.5:8; v=exp(v); v=v*range/max(v);
[X Y]=meshgrid(ms,cs);
contour(X,Y,Z1,v,'LineWidth',2);
colormap([0 0 0]);
set(gca,'Linewidth',2); set(gca,'FontSize',20);
xlabel('Slope, {\it m}'); ylabel('Intercept, {\it c}');
figure(3);
grid on;
set(gca,'LineWidth',2);
surfl(X,Y,Z1)
shading interp
set(gca,'Linewidth',2); set(gca,'FontSize',20);
xlabel('{\it m}');
ylabel('{\it c}');
zlabel('F');
view(-13,30);
colormap(gray);
set(gca,'XLim',[-0.5 1.5],'FontName','Ariel');
set(gca,'YLim',[2 4],'FontName','Ariel');
return