%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Widely Applicable Bayesian Information Criterion (WBIC) %%%
%%% Made by Sumio Watanabe
%%% This program shows a method how to use WBIC
%%% in reduced rank regression.
%%% Sumio Watanabe, ``A widely applicable Bayesian
%%% information crterion," Journal of Machine Learning Research,
%%% Vol.14, (Mar), pp.867-897, 2013.
%%%%%%%%%%%%%%%%%%%% Clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf;
clear;
close all;
%%%%%%%%%%%%%%%%%%%% Random Environment set %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% You can use any natural numbers for RAND_1 and RAND_2
%%% True parameter and training samples are randomly generated.
RAND_1=10; %%% random seed
RAND_2=10; %%% random seed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Constants in Reduced Rank Regression Y=B*A*X+noise %%%%%%%%%%%%%%
MM=6; %%% Input Dimension
HMAX=6; %%% Hidden Dimension
NN=6; %%% Output Dimension
H0=3; %%% True Hidden Dimension
SD1=3; %%% Standard Deviation of Input Distribution
SD2=0.1; %%% Standard Deviation of Output Noise
SD3=10.0; %%% Standard Deviation of Prior Distribution
SD6=0.2; %%% Standard Deviation of True Parameter Making
%%%%%%% Constants in Metropolis Method %%%%%%%%%%%%%%%%%
KK=2000; %%% Number of Parameters from Posterior Distribution
n=500; %%% Number of training samples;
BURNIN=50000; %%% Burn-in in Metroplois method
INTER=100; %%% sampling interval in Metroplois method
MONTEC=BURNIN+KK*INTER; %%% Total Samling in Metroplis method
BETA=1/log(n); %%% inverse temperature
SMALLVAL=0.5; %%% Constant for calculation of RLCT
SD4=0.0012; %%% Standard Deviation of Monte Carlo Step for A
SD5=0.0012; %%% Standard Deviation of Monte Carlo Step for B
%%% SD4 and SD5 are determined so that 0.1< acceptance probability <0.9.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Random seeds are set %%%%%%%%%%%%%%%%%%%%%%%%%%%%
rand('state',RAND_1);
randn('state',RAND_2);
%%%%%%%%%%%%%%%%%%%%%%%%% Ture parameter is determined %%%%%%%%%%%%%%%%
A0=SD6*randn(H0,MM);
B0=SD6*randn(NN,H0);
%%%%%%%%%%%%%%%%%%%%%%%%%% Input and Output are determined %%%%%%%%%%%%%%
X = SD1*randn(MM,n); %%% random inputs
Y = B0*A0*X + SD2*randn(NN,n); %%% random outputs
%%%%%%%%%%%%%%%%%%% Functions for Likelihood and prior %%%%%%%%%%%%%%%%%%%%
LOGL=@(A,B)((1/(2*SD2*SD2)*(trace((Y-B*A*X)*(Y-B*A*X)')))+NN*n*log(SD2));
PRIO=@(A,B)(1/(2*SD3*SD3)*(trace(A*A')+trace(B*B')));
%%% In reduced rank regression,
%%% LOGL = - \sum log p_i
%%% = \sum (1/2sigma^2)|y_i-BAx_i|^2 + NN*n/2*log(2*pi*sigma^2)
%%% = (1/2sigma^2) tr( (Y-BAX)(Y-BAX)' ) + NN*n*log(sigma) + const
%%% PRIO = - log (prior distribution)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LLL=zeros(1,KK); %%%%% Log Likelihood for parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Model selection by WBIC in Reduced Rank regression %g-%g-%g. True rank = %g.\n',MM,HMAX,NN,H0);
fprintf('Number of Training samples = %g. Number of Metropolis Parameters = %g.\n',n,KK);
%%%%%%%%%%% Model Selection Start %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for HH=1:1:HMAX
%%%%%%%%%%%%%%%%%%% Initial parameters %%%%%%%%%%%%%%%%%%%%%%%%%
A=zeros(HH,MM);
B=zeros(NN,HH);
%%%%%%%%%%%%%%%%%%%%% Metropolis preparation %%%%%%%%%%%%%%%%%%%%%%%%%
ENERGY1=LOGL(A,B);
HAMILTON1=BETA*ENERGY1+PRIO(A,B);
rec=0;
acceptance=0;
maxlll=0;
%%%%%%%%%%%%%%%%%%%%%% Metropolis BEGIN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for t=1:1:MONTEC
%%%%%%%%%%%%%%%%%%%%%% Metropolis Step %%%%%%%%%%%%%%%
AA=A+SD4*randn(HH,MM);
BB=B+SD5*randn(NN,HH);
ENERGY2=LOGL(AA,BB);
HAMILTON2=BETA*ENERGY2+PRIO(AA,BB);
DELTA=HAMILTON2-HAMILTON1;
%%%%%%%%%%%%%%%%%%%% Accept or Reject %%%%%%%%%%%%%%%%%%%%
if(exp(-DELTA)>rand(1,1))
A=AA;
B=BB;
HAMILTON1=HAMILTON2;
ENERGY1=ENERGY2;
if(t>BURNIN)
acceptance=acceptance+1;
end
end
%%%%%%%%%%%%%%%%%%%%%% Record Likelihood %%%%%%%%%%%%%%%%%%
if(mod(t,INTER)==0&t>BURNIN)
rec=rec+1;
LLL(rec)=ENERGY1;
if((t==1)|(ENERGY1>maxlll))
maxlll=ENERGY1;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Metropolis END %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Estimate RLCT BEGIN %%%%%%%%%%%%%%%%%%%%%%%%%%%%
sum1 = mean(LLL);
sum2 = mean(exp(log(LLL)-SMALLVAL*BETA*(LLL-maxlll)));
sum3 = mean(exp(-SMALLVAL*BETA*(LLL-maxlll)));
sum2 = sum2/sum3;
lambda2=(sum1-sum2)/((1-(1/(1+SMALLVAL)))*log(n));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Estimate RLCT End %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Theory RCLT of reduced rank regression BEGIN %%%%%%%%%%%%%%%
lambda1=(2*(HH+H0)*(MM+NN)-(MM-NN)*(MM-NN)-(HH+H0)*(HH+H0))/8;
if(mod(MM+NN+HH+H0,2)==1)
lambda1=(2*(HH+H0)*(MM+NN)-(MM-NN)*(MM-NN)-(HH+H0)*(HH+H0)+1)/8;
end
if(HH