forked from burakbayramli/books
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlinregDemo.m
109 lines (87 loc) · 2.69 KB
/
linregDemo.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
seed = 0;
randn('state', seed);
rand('state', seed);
useRBF = 1;
doSave = 0;
longRange = 0;
trainStep = 0.1;
%trainStep = 0.5;
if longRange
xTrainRaw = [-9:trainStep:-8, 3:trainStep:3.5]'+0;
xTestRaw = [-10:0.1:10]'+0;
else
xTrainRaw = [-1:trainStep:-0.5,3:trainStep:3.5]'+0;
xTestRaw = [-2:0.1:4.5]'+0;
end
yTrain = feval(@fun, xTrainRaw) + randn(size(xTrainRaw,1),1)*10;
if useRBF
lo = min( min(xTrainRaw), min(xTestRaw) );
hi = min( max(xTrainRaw), max(xTestRaw) );
RBFcenters = linspace(lo, hi, 5)'
%RBFcenters = linspace(-5, 5, 5)';
xTrain = rbfBasis(xTrainRaw, RBFcenters);
else
xTrain = polyBasis(xTrainRaw, 2);
end
yTestOpt = feval(@fun, xTestRaw);
if useRBF
xTest = rbfBasis(xTestRaw, RBFcenters);
else
xTest = polyBasis(xTestRaw, 2);
end
[xTrain, muTrain, sigmaTrain] = standardizeOnes(xTrain);
w = xTrain\ytrain;
xTest = standardizeOnes;
yPred =
[yPred, yPred_var] = bayes_fwd(mu_w, Sigma_w, beta, xTest);
[yPredTrain, yPredTrain_var] = bayes_fwd(mu_w, Sigma_w, beta, xTrain);
% Plot predictive distribution
figure(1); clf
hold off
subplot(111)
errorbar(xTestRaw(:,1), yPred, sqrt(yPred_var), 'kx-');
hold on
plot(xTestRaw(:,1),yTestOpt,'bx-');
h = errorbar(xTrainRaw(:,1), yPredTrain, sqrt(yPredTrain_var), 'gx');
h = plot(xTrainRaw(:,1), yTrain, 'ro');
set(h, 'linewidth', 3)
grid on
if doSave
fname = sprintf('bayes_lin_regr_demo_frank_rbf%d.eps', useRBF);
folder = 'C:\kmurphy\Teaching\stat406-spring06\Book\figures';
print(gcf, '-depsc', fullfile(folder, fname));
fname = sprintf('bayes_lin_regr_demo_frank_rbf%d.jpg', useRBF);
folder = 'C:\kmurphy\Teaching\stat406-spring06\Book\figures';
print(gcf, '-djpeg', fullfile(folder, fname));
end
% Plot samples from the posterior
figure(2); clf
plot(xTestRaw(:,1),yTestOpt,'bx-'); hold on
nsamples = 10;
for s=1:nsamples
w = mvnrnd(mu_w, Sigma_w, 1);
Y_mean = xTest * w(:);
plot(xTestRaw(:,1), Y_mean, 'kx-');
end
if doSave
fname = sprintf('bayes_lin_regr_demo_frank_samples_rbf%d.eps', useRBF);
folder = 'C:\kmurphy\Teaching\stat406-spring06\Book\figures';
print(gcf, '-depsc', fullfile(folder, fname));
fname = sprintf('bayes_lin_regr_demo_frank_samples_rbf%d.jpg', useRBF);
folder = 'C:\kmurphy\Teaching\stat406-spring06\Book\figures';
print(gcf, '-djpeg', fullfile(folder, fname));
end
%%%%%%%%
function f = fun(x) %target function
f = x.^4;
%%%%%%%%
function [mu_w, Sigma_w] = bayesian_update(mu_w, Sigma_w, beta, X, Y)
Sigma_old_inv = inv(Sigma_w);
Sigma_w = pinv(Sigma_old_inv + (beta * X' * X));
mu_w = (Sigma_w*Sigma_old_inv) * mu_w + (beta * Sigma_w * X' * Y);
%%%%%%%%
function [Y_mean, Y_var] = bayes_fwd(mu_w, Sigma_w, beta, X)
Y_mean = X * mu_w;
for i=1:size(X,1)
Y_var(i,1) = 1/beta + X(i,:) * Sigma_w * X(i,:)';
end