forked from ExaScience/bpmf
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathbayespmf.m
162 lines (129 loc) · 4.63 KB
/
bayespmf.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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
% Version 1.000
%
% Code provided by Ruslan Salakhutdinov
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, express or
% implied. As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application. All use of these programs is entirely at the user's own risk.
rand('state',0);
randn('state',0);
if restart==1
restart=0;
epoch=1;
maxepoch=50;
iter=0;
num_m = 3952;
num_p = 6040;
num_feat = 10;
% Initialize hierarchical priors
beta=2; % observation noise (precision)
mu_u = zeros(num_feat,1);
mu_m = zeros(num_feat,1);
alpha_u = eye(num_feat);
alpha_m = eye(num_feat);
% parameters of Inv-Whishart distribution (see paper for details)
WI_u = eye(num_feat);
b0_u = 2;
df_u = num_feat;
mu0_u = zeros(num_feat,1);
WI_m = eye(num_feat);
b0_m = 2;
df_m = num_feat;
mu0_m = zeros(num_feat,1);
load moviedata
mean_rating = mean(train_vec(:,3));
ratings_test = double(probe_vec(:,3));
pairs_tr = length(train_vec);
pairs_pr = length(probe_vec);
fprintf(1,'Initializing Bayesian PMF using MAP solution found by PMF \n');
makematrix
load pmf_weight
err_test = cell(maxepoch,1);
w1_P1_sample = w1_P1;
w1_M1_sample = w1_M1;
clear w1_P1 w1_M1;
% Initialization using MAP solution found by PMF.
%% Do simple fit
mu_u = mean(w1_P1_sample)';
d=num_feat;
alpha_u = inv(cov(w1_P1_sample));
mu_m = mean(w1_M1_sample)';
alpha_m = inv(cov(w1_P1_sample));
count=count';
probe_rat_all = pred(w1_M1_sample,w1_P1_sample,probe_vec,mean_rating);
counter_prob=1;
end
for epoch = epoch:maxepoch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Sample from movie hyperparams (see paper for details)
N = size(w1_M1_sample,1);
x_bar = mean(w1_M1_sample)';
S_bar = cov(w1_M1_sample);
WI_post = inv(inv(WI_m) + N/1*S_bar + ...
N*b0_m*(mu0_m - x_bar)*(mu0_m - x_bar)'/(1*(b0_m+N)));
WI_post = (WI_post + WI_post')/2;
df_mpost = df_m+N;
alpha_m = wishrnd(WI_post,df_mpost);
mu_temp = (b0_m*mu0_m + N*x_bar)/(b0_m+N);
lam = chol( inv((b0_m+N)*alpha_m) ); lam=lam';
mu_m = lam*randn(num_feat,1)+mu_temp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Sample from user hyperparams
N = size(w1_P1_sample,1);
x_bar = mean(w1_P1_sample)';
S_bar = cov(w1_P1_sample);
WI_post = inv(inv(WI_u) + N/1*S_bar + ...
N*b0_u*(mu0_u - x_bar)*(mu0_u - x_bar)'/(1*(b0_u+N)));
WI_post = (WI_post + WI_post')/2;
df_mpost = df_u+N;
alpha_u = wishrnd(WI_post,df_mpost);
mu_temp = (b0_u*mu0_u + N*x_bar)/(b0_u+N);
lam = chol( inv((b0_u+N)*alpha_u) ); lam=lam';
mu_u = lam*randn(num_feat,1)+mu_temp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start doing Gibbs updates over user and
% movie feature vectors given hyperparams.
for gibbs=1:2
fprintf(1,'\t\t Gibbs sampling %d \r', gibbs);
%%% Infer posterior distribution over all movie feature vectors
count=count';
for mm=1:num_m
fprintf(1,'movie =%d\r',mm);
ff = find(count(:,mm)>0);
MM = w1_P1_sample(ff,:);
rr = count(ff,mm)-mean_rating;
covar = inv((alpha_m+beta*MM'*MM));
mean_m = covar * (beta*MM'*rr+alpha_m*mu_m);
lam = chol(covar); lam=lam';
w1_M1_sample(mm,:) = lam*randn(num_feat,1)+mean_m;
end
%%% Infer posterior distribution over all user feature vectors
count=count';
for uu=1:num_p
fprintf(1,'user =%d\r',uu);
ff = find(count(:,uu)>0);
MM = w1_M1_sample(ff,:);
rr = count(ff,uu)-mean_rating;
covar = inv((alpha_u+beta*MM'*MM));
mean_u = covar * (beta*MM'*rr+alpha_u*mu_u);
lam = chol(covar); lam=lam';
w1_P1_sample(uu,:) = lam*randn(num_feat,1)+mean_u;
end
end
probe_rat = pred(w1_M1_sample,w1_P1_sample,probe_vec,mean_rating);
probe_rat_all = (counter_prob*probe_rat_all + probe_rat)/(counter_prob+1);
counter_prob=counter_prob+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Make predictions on the validation data %%%%%%%
temp = (ratings_test - probe_rat_all).^2;
err = sqrt( sum(temp)/pairs_pr);
iter=iter+1;
overall_err(iter)=err;
fprintf(1, '\nEpoch %d \t Average Test RMSE %6.4f \n', epoch, err);
end