anneal.m

function params = annealit(initparams, type)
%annealit:  a simple annealing routine

%Annealing is based on a simple idea of jumping around a local parameter space looking for the minima
%A proposed parameter set is accepted with probability exp(-(energy_old_params - energy_new_params)/temperature) the 'energy' (in this example energy = neg log likelihood)
%  The effect is that if energy_new is lower then the new param set is accepted but if
%  energy_new is higher it accepted with a probability that depends on how much higher it is
%  Temperature affects how this probability depends on the energy difference.
%  The idea is to slowly lower the temperature to get to the lowest energy: voila, simulated annealing
%  Surprisingly this simple little idea often works

%initparams is composed of initparams.free and initparams.fixed; it's a vector
%fittomodel is some model that is called with params and returns -log(L(data|params))

%iter is the number of annealing calls to make
if(nargin == 1) type = 2; end
iter=10000;
limunif = .1; %10% up or down
%T is the temperature; a somewhat vague concept with basically indicates how much the annealing jumps around
T = 10;

%initialize things
energy_old=fittomodel(initparams); 
param=initparams;
srecord = zeros(iter,1); %just a holding bin

%run iter iterations of the annealing routine
for(j = 1:iter)
   %in this routine, I go through the free params one by one
   %one can randomly pick a completely new set, but then it is quite unlikely that one
   %will choose a new set with lower energy
   for(i = 1:length(params.free)) 
      proposed_params=params;
      proposed_params.free=params.free;
      switch type %different ways you can choose new params; I like #2
         case 1
            proposed_params.free(i)=unifrnd( (1-limunif)*params.free(i),(1+limunif)*params.free(i) );
         case 2
            proposed_params.free(i)=unifrnd( -limunif,limunif )+params.free(i);
         otherwise
            display('type must be 1 or 2'); break
      end
      energy_new=fittomodel(proposed_params);
      %switch probability
      swprob = unifrnd(0,1);
      ener = exp(-(energy_new - energy_old)/T);
      if(swprob < ener) %accept
         params.free = proposed_params.free;
         energy_old = energy_new;
      end
   end
   srecord(j) = fittomodel(params);
   disp(srecord(j));
   if(mod(j,100)==0) 
      display(['iter = ' num2str(j) '   T = ' num2str(T)]);
      display('Now lower T slowly to converge to global minima');
      display('Yes, manually do this by typing T = y; y being something lower than Ts current value');
      keyboard; 
   end
end
plot(srecord) %this will help you see if the annealing is working