-
Notifications
You must be signed in to change notification settings - Fork 33
/
gibbs.m
123 lines (79 loc) · 2.39 KB
/
gibbs.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
119
120
121
122
123
% ---------------------------
% Simulation of a Gibbs sampling
% written by Richard Xu
% March 2013
% --------------------------
% clear;
% m = 4; k=3; n = 5; A = randn(m,k); B = randn(k,n); % values
% a1 = A * B; % method 1
% display(a1);
%
% a2 = zeros(m,n);
% for i=1:k; % method 2
% a2 = a2+ A(:,i) * B(i,:);
% end;
% display(a2);
%
% return;
clc;
clf;
clear;
mu = [2 3];
Sigma = [ 3 2; 2 5];
figure(1);
N = 2000;
init_sample = [ 1 1];
directX = mvnrnd(mu,Sigma,N);
minX = min(directX(:,1));
minY = min(directX(:,2));
maxX = max(directX(:,1));
maxY = max(directX(:,2));
figure(1);
subplot(2,2,1);
plot(directX(:,1),directX(:,2),'.','color',[1 0 0]);
axis([minX maxX minY maxY]);
subplot(2,2,2);
autocorr(mvnpdf(directX,mu,Sigma));
HIDE_GIBBS_STEP = false;
if HIDE_GIBBS_STEP == true
X = zeros(N,2);
X(1,:) = init_sample;
for i = 2:N
mu_1 = mu(1) + Sigma(1,2)/Sigma(2,2)*( X(i-1,2) - mu(2));
Sigma_1 = Sigma(1,1) - Sigma(1,2)/Sigma(2,2)*Sigma(2,1);
X(i,1) = randn * sqrt(Sigma_1) + mu_1;
mu_2 = mu(2) + Sigma(2,1)/Sigma(1,1)*(X(i,1) - mu(1));
Sigma_2 = Sigma(2,2) - Sigma(2,1)/Sigma(1,1)*Sigma(1,2);
X(i,2) = randn * sqrt(Sigma_2) + mu_2;
end
subplot(2,2,3);
plot(X(:,1),X(:,2),'.');
axis([minX maxX minY maxY]);
%---------------------------------------
else
subplot(2,2,3);
X = zeros(N,2);
X(1,:) = [1 1];
for i = 2:N
mu_1 = mu(1) + Sigma(1,2)/Sigma(2,2)*( X(i-1,2) - mu(2));
Sigma_1 = Sigma(1,1) - Sigma(1,2)/Sigma(2,2)*Sigma(2,1);
X(i,1) = randn * sqrt(Sigma_1) + mu_1;
plot( [X(i-1,1) X(i,1)], [X(i-1,2) X(i-1,2)],'LineWidth',1);
axis([minX maxX minY maxY]);
hold on;
waitforbuttonpress;
%waitfor(100);
mu_2 = mu(2) + Sigma(2,1)/Sigma(1,1)*(X(i,1) - mu(1));
Sigma_2 = Sigma(2,2) - Sigma(2,1)/Sigma(1,1)*Sigma(1,2);
X(i,2) = randn * sqrt(Sigma_2) + mu_2;
plot( [X(i,1) X(i,1)], [X(i-1,2) X(i,2)],'LineWidth',1);
axis([minX maxX minY maxY]);
hold on;
waitforbuttonpress;
%waitfor(100);
end
plot(X(:,1),X(:,2),'.');
end
subplot(2,2,4);
autocorr(mvnpdf(X,mu,Sigma));