0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 2 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 3 % Metropolis algorithm
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 4 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 5 % M Nofrarias 15-06-09
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 6 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 7 % $Id: mhsample_td.m,v 1.1 2010/11/26 13:48:48 luigi Exp $
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 8 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 10
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 11
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 12 function smpl = mhsample_td(model,in,out,cov,number,limit,parnames,Tc,xi,xo,search,jumps,parplot,dbg_info,inNames,outNames,Noise,cutbefore,cutafter)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 13
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 14 import utils.const.*
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 15
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 16 % compute mean and range
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 17 mn = mean(limit);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 18 range = diff(limit);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 19 % initialize
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 20 acc = [];
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 21 nacc = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 22 nrej = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 23 loop = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 24 oldacc = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 25 oldrej = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 26 chgsc = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 27 nparnames = length(parnames);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 28
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 29
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 30 % switch depending on model class (initial likelihood)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 31 switch class(model)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 32 case 'matrix'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 33 % loglk1 = utils.math.loglikehood_matrix(xo,in,out,nse,model,parnames);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 34 % loglk1 = utils.math.loglikehood_matrix(xo,in,out,nse,model,parnames,inModel,outModel);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 35 case 'smodel'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 36 % loglk1 = utils.math.loglikehood(xo,in,out,nse,model,parnames);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 37 case 'ssm'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 38 loglk1 = utils.math.loglikehood_ssm_td(xo,in,out,parnames,model,inNames,outNames,Noise,'cutbefore',cutbefore,'cutafter',cutafter);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 39 otherwise
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 40 error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 41 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 42
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 43 % accept the first sample if no search active
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 44 if ~search
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 45 smpl(1,:) = xo;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 46 nacc = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 47 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 48 smpl = [];
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 49 nacc = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 50 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 51 nrej = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 52
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 53 utils.helper.msg(msg.IMPORTANT, 'Starting Monte Carlo sampling', mfilename('class'), mfilename);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 54
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 55
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 56 hjump = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 57 % main loop
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 58 while(nacc<number)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 59 % jumping criterion during search phase
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 60 % - 2-sigma by default
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 61 % - 10-sigma jumps at mod(10) samples
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 62 % - 100-sigma jumps at mod(25) samples
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 63 % - 1000-sigma jumps at mod(100) samples
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 64 if search
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 65 if nacc <= Tc(1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 66 if(mod(nacc,10) == 0 && mod(nacc,25) ~= 0 && mod(nacc,100) ~= 0 && hjump ~= 1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 67 hjump = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 68 xn = mvnrnd(xo,jumps(2)^2*cov);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 69 elseif(mod(nacc,25) == 0 && mod(nacc,100) ~= 0 && hjump ~= 1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 70 hjump = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 71 xn = mvnrnd(xo,jumps(3)^2*cov);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 72 elseif(mod(nacc,100) == 0 && hjump ~= 1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 73 hjump = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 74 xn = mvnrnd(xo,jumps(4)^2*cov);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 75 % xn = xo + range/2+rand(size(xo));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 76 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 77 hjump = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 78 xn = mvnrnd(xo,jumps(1)^2*cov);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 79 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 80 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 81 xn = mvnrnd(xo,cov);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 82 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 83 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 84 xn = mvnrnd(xo,cov);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 85 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 86
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 87 % check if out of limits
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 88 if( any(xn < limit(1,:)) || any(xn > limit(2,:)))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 89 logprior = inf;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 90 loglk2 = inf;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 91 betta = inf;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 92 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 93 % compute annealing
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 94 if ~isempty(Tc)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 95 if nacc <= Tc(1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 96 betta = 1/2 * 10^(-xi*(1-Tc(1)/Tc(2)));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 97 elseif Tc(1) < nacc && nacc <= Tc(2)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 98 betta = 1/2 * 10^(-xi*(1-nacc/Tc(2)));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 99 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 100 betta = 1/2;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 101 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 102 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 103 betta = 1/2;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 104 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 105 % switch depending on model class (new likelihood)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 106 switch class(model)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 107 case 'matrix'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 108 % loglk2 = utils.math.loglikehood_matrix(xn,in,out,nse,model,parnames,inModel,outModel);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 109 case 'smodel'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 110 % loglk2 = utils.math.loglikehood(xn,in,out,nse,model,parnames);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 111 case 'ssm'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 112 loglk2 = utils.math.loglikehood_ssm_td(xn,in,out,parnames,model,inNames,outNames,Noise,'cutbefore',cutbefore,'cutafter',cutafter);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 113 otherwise
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 114 error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 115 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 116 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 117
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 118 % likelihood ratio
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 119 logr = betta*(loglk2 - loglk1) ;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 120 % decide if sample is accepted or not
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 121 if logr < 0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 122 xo = xn;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 123 nacc = nacc+1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 124 smpl(nacc,:) = xn;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 125 loglk1 =loglk2;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 126 if dbg_info
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 127 utils.helper.msg(msg.IMPORTANT, sprintf('acc.\t loglik: %d -> %d betta: %d ratio: %d',loglk1,loglk2,betta,logr), mfilename('class'), mfilename);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 128 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 129 elseif rand(1) > (1 - exp(-logr))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 130 xo = xn;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 131 nacc = nacc+1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 132 smpl(nacc,:) = xn;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 133 loglk1 =loglk2;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 134 if dbg_info
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 135 utils.helper.msg(msg.IMPORTANT, sprintf('acc.\t loglik: %d -> %d betta: %d ratio: %d',loglk1,loglk2,betta,logr), mfilename('class'), mfilename);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 136 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 137 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 138 nrej = nrej+1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 139 if dbg_info
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 140 utils.helper.msg(msg.IMPORTANT, sprintf('rej.\t loglik: %d -> %d betta: %d ratio: %d',loglk1,loglk2,betta,logr), mfilename('class'), mfilename);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 141 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 142 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 143
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 144 % display and save
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 145 if(mod(nacc,100) == 0 && nacc ~= oldacc)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 146 updacc = nacc-oldacc;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 147 updrej = nrej-oldrej;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 148 ratio(nacc/10,:) = updacc/(updacc+updrej);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 149 utils.helper.msg(msg.IMPORTANT, sprintf('accepted: %d rejected : %d acc. rate: %4.2f',nacc,updrej,ratio(end)), mfilename('class'), mfilename);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 150 for i = 1:numel(parnames)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 151 fprintf('### Parameters: %s = %d \n',parnames{i},xn(i))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 152 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 153 oldacc = nacc;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 154 oldrej = nrej;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 155 save('chain.txt','smpl','-ASCII')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 156 save('acceptance.txt','ratio','-ASCII')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 157 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 158
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 159 % plot
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 160 if ((mod(nacc,100) == 0) && (nacc ~= 0) && ~isempty(parplot))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 161 for i = 1:numel(parplot)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 162 figure(parplot(i))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 163 plot(smpl(:,parplot(i)))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 164 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 165 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 166
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 167 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 168 end