-
Notifications
You must be signed in to change notification settings - Fork 0
/
exercise6_L4_monte.m
118 lines (83 loc) · 2.91 KB
/
exercise6_L4_monte.m
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
111
112
113
114
115
116
117
118
figure(1); clf; figure(2); clf
x=[2 2.5 2.5 2.75 3 3 3];
y=[89 97 91 98 100 104 97];
X=[ones(size(x')) x'];
y=y';
b=inv(X'*X)*X'*y
P=2; % two parameters
N=length(x); %number of observatrions
nu=N-P; % degrees of freedom. no. obs-no.parameters
bestmodel=X*b;
s2=sum((bestmodel-y).^2)./nu; s=sqrt(s2)
tvalue=tinv(0.975,nu);
seB=s*sqrt(diag(inv(X'*X)))*tvalue
w=0:0.1:2.5*pi; % for plotting the elippse
[Q,R]=qr(X); R1 = R(1:2,1:2); invR1=inv(R1); Fvalue=finv(0.95,P,N-P);
for i=1:length(w)
scalar=sqrt(P*s2*Fvalue);
BETA(:,i)=b+scalar*invR1*([cos(w(i)); sin(w(i))]);
end
%plot(BETA(1,:),BETA(2,:),'linewidth',2)
%set(gca,'fontsize',11,'linewidth',2)
%xlabel('intercept'); ylabel('slope')
% mc test
residuals=y-bestmodel; sestimate=std(residuals);
% need to make simulated dataset
count = 1:length(y); % an integer vector with an entry in order per each point in the data
% replacing all data seems to work better
%noreplace=round(length(y)/3); % replace 1/3 of the datapoints
noreplace=round(length(y)/1); % replace all of the datapoints
indexreplace = randsample(count,noreplace); % pick at random wich of the 1/3 of the points to replace.
Ynew=y; % initialize the simulated data just as the original data.
for i=1:length(indexreplace)
Ynew(indexreplace(i))=bestmodel(indexreplace(i))+randn(1,1)*sestimate;
end
%plot(x,y,'ko',x,Ynew,'k.')
% full MC
for j=1:1000
Ynew=y; % initialize the simulated data just as the original data.
for i=1:length(indexreplace)
Ynew(indexreplace(i))=bestmodel(indexreplace(i))+randn(1,1)*sestimate;
end
MCbetas(:,j)=inv(X'*X)*X'*Ynew;
fit(:,j)=X*MCbetas(:,j);
end
MCslope=MCbetas(2,:); MCintercepts=MCbetas(1,:);
mean(MCslope)
std(MCslope)
figure(1)
subplot(211);
[counts,centres]=hist(MCslope,10);
hist(MCslope,10)
hold on; plot([b(2) b(2)],[0 max(counts)],'k','linewidth',2)
plot([b(2)-seB(2) b(2)-seB(2)],[0 max(counts)],'k--','linewidth',2)
plot([b(2)+seB(2) b(2)+seB(2)],[0 max(counts)],'k--','linewidth',2)
axis([b(2)-seB(2)*1.2 b(2)+seB(2)*1.2 0 max(counts)*1.1])
subplot(212)
[counts,centres]=hist(MCintercepts,10);
hist(MCintercepts,10)
hold on; plot([b(1) b(1)],[0 max(counts)],'k','linewidth',2)
plot([b(1)-seB(1) b(1)-seB(1)],[0 max(counts)],'k--','linewidth',2)
plot([b(1)+seB(1) b(1)+seB(1)],[0 max(counts)],'k--','linewidth',2)
axis([b(1)-seB(1)*1.2 b(1)+seB(1)*1.2 0 max(counts)*1.1])
% plot data and swarm of lines
figure(2)
plot(x,bestmodel,x,y,'ko')
hold on
plot(x,fit)
% add the dashed lines
Fvalue=finv(0.95,P,N-P);
xplot=1.8:0.1:3.2; xplot=xplot';
XPLOT=[ones(size(xplot)) xplot];
model=XPLOT*b;
plot(x,y,'ko','markersize',4,'markerfacecolor','b')
set(gca,'linewidth',2,'fontsize',11)
hold on
plot(xplot,model,'k','linewidth',2)
for i=1:length(xplot)
xh=[1; xplot(i)];
upper(i)=xh'*b+s*sqrt(xh'*inv(X'*X)*xh)*sqrt(P*Fvalue);
lower(i)=xh'*b-s*sqrt(xh'*inv(X'*X)*xh)*sqrt(P*Fvalue);
end
plot(xplot,upper,'k--')
plot(xplot,lower,'k--')