MCMC Optimization

Example  1: continuous function with a single maximum

Here is a function to maximize:

f[x_]:=Exp[-(x^2)]
Plot[f[x],{x,-3,3},PlotRange->All]

Here is the Markov chain Monte Carlo code. This module requires the following parameters:

δ: a nonnegative radius used to select proposal transitions.

β: a nonnegative weight parameter

length: the length of Markov chain to generate

start: an initial value for the Markov chain

mcmc[δ_,β_,length_,start_,f_]:=Module[
(* local variables *)
{chain=Table[0,length], (* list for storing the Markov chain *)
currX=start,propX,Δf,ρ
},
(* store the initial value in the chain *)
chain[[1]]=currX;

(* generate the Markov chain *)
Do[
(* select a proposal transition *)
currX=chain[[m-1]];
propX=;

(* compute the acceptance ratio *) 
Δf=f[propX]-f[currX];
ρ=;

(* decide whether to accept the proposed transition *) 
If[RandomReal[]<ρ,
currX=propX
];

(* store the next state *) 
chain[[m]]=currX
,{m,2,length}];

(* return the chain *)
chain
]

Plot  the Markov chain:

chain=mcmc[]
ListLinePlot[chain,PlotRange->All]

Questions

Generate Markov chains for various choices of δ and β. How do these parameters affect the Markov chain?

What proportion of the time is the proposal transition accepted? How does this depend on δ and β?

How quickly does the Markov chain converge to 0? How does this depend on δ and β?

Example 2: continuous function with three local maxima

Now consider the following function.

g[x_]:=0.25(Exp[-10(((x-1.5))^2)]+Exp[-10(((x+1.5))^2)])+Exp[-10(x^2)]
Plot[g[x],{x,-3,3},PlotRange->All]

Use Markov chain Monte Carlo to search for the global maximum of g(x).

Questions

Generate and plot various Markov chains. How would you describe the Markov chains? What behaviors do you notice?

How do δ and β affect the proportion of times that the proposal transition is accepted?

Does the Markov chain reliably identify the global maximum of g(x)? Why or why not?

Can you more reliably identify the global maximum of g(x) by starting with a small β and gradually increasing β?

Discrete State Space: Integer Partitions

A partition of an integer N is a way of writing N as a sum of positive integers. For example, N=4 has five distinct partitions:

4
3+1
2+2
2+1+1
1+1+1+1

Does there exist a partition of N=1000 consisting of some of the following 50 integers?

target=1000;
integers={26,52,44,44,76,30,46,34,27,90,97,91,43,14,36,86,93,90,11,14,62,71,65,47,72,13,79,69,95,10,50,37,69,27,13,27,78,50,53,12,94,95,81,22,36,74,54,28,18,64};

Use MCMC methods to try to find such a partition of N=1000.

First, we have to decide how to represent a state:

currState=

We also need a function to maximize:

f[state_]:=

Now we can run the simulation:

(* initialize values *)
maxIter=100;


(* main loop *)
Do[
(* select a proposed next state *)


(* select the next β in the annealing schedule *)


(* compute the acceptance ratio *) 
Δf=f[propState]-f[currState];
ρ=;

(* decide whether to accept the proposed transition *) 
If[RandomReal[]<ρ,
currState=propState
];

(* if a partition has been found, then quit *)

,maxIter]

Questions

Can you find a subset of the given integers that adds up to N=1000? Is this subset unique?

How reliably can you find such a subset?

What annealing schedule works best?

Try other targets and other sets of integers. Can you replicate your results consistently?