rjMCMCsa.rar_As One_Bayesian_de Freitas_neural matlab _rjMCMCsa

function [k,mu,alpha,ypred,ypredv,post] = rjsann(x,y,chainLength,Ndata,bFunction,par,xv,yv); % PURPOSE : Computes the parameters and number of parameters of a radial basis function (RBF) % network using the reversible jump MCMC simulated annealing algorithm. Please have a % look at the paper first. % INPUTS : - x : Input data. % - y : Target data. % - chainLength: Number of iterations of the Markov chain. % - Ndata: Number of time steps in the training data set. % - bFunction: Type of basis function. % - par: Record of simulation parameters (see defaults). % - {xv,yv}: Validation data (optional). % AUTHOR : Nando de Freitas - Thanks for the acknowledgement :-) % DATE : 21-01-99 % CHECK INPUTS AND SET DEFAULTS: % ============================= if nargin < 5, error('Not enough input arguments.'); end if ((nargin==5) | (nargin==7)), if nargin == 5, Validation = 0; else Validation = 1; end; hyper.a = 2; % Hyperparameter for delta. hyper.b = 10; % Hyperparameter for delta. hyper.e1 = 0.0001; % Hyperparameter for nabla. hyper.e2 = 0.0001; % Hyperparameter for nabla. hyper.v = 0; % Hyperparameter for sigma hyper.gamma = 0; % Hyperparameter for sigma. kMax = 50; % Maximum number of basis. arbC = 0.25; % Constant for birth and death moves. doPlot = 1; % To plot or not to plot? Thats ... sigStar = .1; % Merge-split parameter. sWalk = .01; Lambda = .5; walkPer = 0.1; Ta = 1; % Cooling parameter (~Initial temperature). Tf = 0.1; % Cooling parameter (Final temperature). elseif ((nargin==6) | (nargin==8)), if nargin == 6, Validation = 0; else Validation = 1; end; Validation == 0 hyper.a = par.a; hyper.b = par.b; hyper.e1 = par.e1; hyper.e2 = par.e2; hyper.v = par.v; hyper.gamma = par.gamma; kMax = par.kMax; arbC = par.arbC; doPlot = par.doPlot; sigStar = par.merge; sWalk = par.sRW; Lambda = par.Lambda; walkPer = par.walkPer; Ta = par.Ta; Tf = par.Tf; else error('Wrong Number of input arguments.'); end; if Validation, [Nv,dv] = size(xv); % Nv = number of test data, dv = dimension of xv. end; [N,d] = size(x) % N = number of data, d = dimension of x. [N,c] = size(y) % c = dimension of y, i.e. number of outputs. if Ndata ~= N, error('input must me N by d and output N by c.'); end % HYPER-PARAMETERS AND INITIALISATION: % =================================== Tb = (Tf-Ta)/chainLength; % Cooling parameter (Slope). Temp = ones(chainLength,1); % SA temperature. post = ones(chainLength,1); % p(centres,k|y). if Validation, ypredv = zeros(Nv,c,chainLength); % Output fit (test set). end; ypred = zeros(N,c,chainLength); % Output fit (train set). k = ones(chainLength,1); % Model order - number of basis. mu = cell(chainLength,1); % Radial basis centres. alpha = cell(chainLength,c); % Radial basis coefficients. if par.criterion == 1, % Model selection criterion: criterion= c+1; % AIC = c+1 elseif par.criterion == 2, criterion= ((c+1)/2) * log(N); % MDL = ((c+1)/2) * log(N) end; % DEFINE WALK INTERVAL FOR MU: % =========================== walk = walkPer*(max(x)-min(x)); walkInt=zeros(d,1); for i=1:d, walkInt(i,1) = (max(x(:,i))-min(x(:,i))) + 2*walk(i); end; % INITIAL CONDITION: (for comparison purposes in exp 1. Delete it otherwise.) % ================= k(1) = 20; % DRAW THE INITIAL RADIAL CENTRES WITHOUT REPLACEMENT: % =================================================== mu{1}=zeros(k(1),d); for i=1:d, mu{1}(:,i)= (min(x(:,i))-walk(i))*ones(k(1),1) + ((max(x(:,i))+walk(i))-(min(x(:,i))-walk(i)))*rand(k(1),1); end % FILL THE REGRESSION MATRIX: % ========================== M=zeros(N,k(1)+d+1); M(:,1) = ones(N,1); M(:,2:d+1) = x; for j=d+2:k(1)+d+1, M(:,j) = feval(bFunction,mu{1}(j-d-1,:),x); end; % COMPUTE THE PREDICTION AT t=0: % ============================= H=zeros(k(1)+1+d,k(1)+1+d,c); F=zeros(k(1)+1+d,c); P=zeros(N,N,c); for i=1:c, H(:,:,i) = inv(M'*M); F(:,i) = H(:,:,i)*M'*y(:,i); P(:,:,i) = eye(N) - M*H(:,:,i)*M'; alpha{1,i} = F(:,i); end; for i=1:c, ypred(:,i,1) = M*alpha{1,i}; end; if Validation, Mv=zeros(Nv,k(1)+d+1); Mv(:,1) = ones(Nv,1); Mv(:,2:d+1) = xv; for j=d+2:k(1)+d+1, Mv(:,j) = feval(bFunction,mu{1}(j-d-1,:),xv); end; for i=1:c, ypredv(:,i,1) = Mv*alpha{1,i}; end; end; % INITIALISE COUNTERS: % =================== Nbirth=0; Ndeath=0; Nupdate=0; aUpdate=0; rUpdate=0; aBirth=0; rBirth=0; aDeath=0; rDeath=0; aMerge=0; rMerge=0; aSplit=0; rSplit=0; aRW=0; rRW=0; match=0; if doPlot, figure(3) clf; end; % ITERATE THE MARKOV CHAIN: % ======================== for t=1:chainLength-1, t=t % thek=k(t) % thepost=post(t) % UPDATE DENOMINATOR FOR NON-HOMOGENEOUS STEP: % ============================================ invH=zeros(k(t)+1+d,k(t)+1+d,c); P=zeros(N,N,c); for i=1:c, invH(:,:,i) = M'*M; P(:,:,i) = eye(N) - M*inv(invH(:,:,i))*M'; end; oldM=M; small = 0; DenRatio = exp(-criterion*k(t)) * (y(:,1)'*P(:,:,1)*y(:,1)+small)^(N/2); for i=2:c, DenRatio = DenRatio * (y(:,i)'*P(:,:,i)*y(:,i)+small)^(N/2); end; % COMPUTE THE CENTRES AND DIMENSION WITH METROPOLIS, BIRTH AND DEATH MOVES: % ======================================================================== decision=rand(1); birth=arbC; death=arbC; if ((decision <= birth) & (k(t)<kMax)), [k,mu,M,match,aBirth,rBirth] = sarBirth(match,aBirth,rBirth,k,mu,M,x,y,t,criterion,bFunction,walkInt,walk); elseif ((decision <= birth+death) & (k(t)>0)), [k,mu,M,aDeath,rDeath] = sarDeath(aDeath,rDeath,k,mu,M,x,y,t,criterion,walkInt); elseif ((decision <= 2*birth+death) & (k(t)<kMax) & (k(t)>1)), [k,mu,M,aSplit,rSplit] = sarSplit(aSplit,rSplit,k,mu,M,x,y,t,bFunction,criterion,sigStar,walk); elseif ((decision <= 2*birth+2*death) & (k(t)>1)), [k,mu,M,aMerge,rMerge] = sarMerge(aMerge,rMerge,k,mu,M,x,y,t,bFunction,criterion,sigStar); else uLambda = rand(1); if ((uLambda>Lambda) & (k(t)>0)) [k,mu,M,match,aRW,rRW] = sarRW(match,aRW,rRW,k,mu,M,x,y,t,criterion,bFunction,sWalk,walk); else [k,mu,M,match,aUpdate,rUpdate] = sarUpdate(match,aUpdate,rUpdate,k,mu,M,x,y,t,criterion,bFunction,walkInt,walk); end; end; % UPDATE NUMERATOR FOR NON-HOMOGENEOUS STEP: % ========================================= invH=zeros(k(t+1)+1+d,k(t+1)+1+d,c); P=zeros(N,N,c); for i=1:c, invH(:,:,i) = M'*M; P(:,:,i) = eye(N) - M*inv(invH(:,:,i))*M'; end; small = 0; NumRatio = exp(-criterion*k(t+1)) * (y(:,1)'*P(:,:,1)*y(:,1)+small)^(N/2); for i=2:c, NumRatio = NumRatio * (y(:,i)'*P(:,:,i)*y(:,i)+small)^(N/2); end; % ANNEALING BIT: % ============= Temp(t) = Ta+Tb*t; % Cooling schedule. U=rand(1); ratio = NumRatio/DenRatio; ratio = ratio^(inv(Temp(t))-1); % Anneal ratio. acceptance = min(1,ratio); if (U<acceptance), mu{t+1} = mu{t+1}; k(t+1) = k(t+1); M=M; else mu{t+1} = mu{t}; k(t+1) = k(t); M=oldM; end; % UPDATE OTHER PARAMETERS WITH GIBBS: % ================================== H=zeros(k(t+1)+1+d,k(t+1)+1+d,c); F=zeros(k(t+1)+1+d,c); P=zeros(N,N,c); for i=1:c, H(:,:,i) = inv(M'*M); F(:,i) = H(:,:,i)*M'*y(:,i); P(:,:,i) = eye(N) - M*H(:,:,i)*M'; alpha{t+1,i} = F(:,i); end; % COMPUTE THE POSTERIOR FOR MONITORING: % ==================================== small = 0; % To avoid numerical problems. posterior