AMOEBA implementation

Here's my code. IDL users, feel free to use this as-is. MATLAB users, feel free to translate it. The point of this assignment is to get the Bayesian periodogram working, not to fiddle with AMOEBA. However, it's a handy tool to have in your belt, so I encourage you to make sure you understand this example of its use.

;===================SINEFIT.PRO======================


function sineorbit, x, per, par
return, par[0] * sin(2*!pi/per * x + par[1]) + par[2]
end

function sinefunc, par
common sinestruct, struct
ymod = sineorbit(struct.x, struct.per, par)
;;; MaxL is just -chi^2 / 2 without the constant. Just use chi as
;;; statistic to be minimized
chi = total((struct.y-ymod)^2/(2.*struct.err)^2)
return, chi
end

function sinefit, x, y, err, period
common sinestruct, struct
struct = {x:x, y:y, err:err, per:period} ;;; create "control" structure

p0 = [(max(y)-min(y))/2., 0., median(y)] ;;; Initital guesses based on data
;;; The scale is the fractional moves in each parameter that AMOEBA will
;;; make while oozing around in maxL (chi^2) space.
scale = [p0[0]*0.05, !pi*0.1, p0[0]*0.05]
;;; Call amoeba and stop after 5% = 0.05 decrease in chi^2 is reached
pbest = amoeba(0.05, function_name='sinefunc', p0=p0, scale=scale)
;;; Calculate best-fitting orbit
ymod = sineorbit(x, period, pbest)
chi = total((struct.y-ymod)^2/(2.*struct.err)^2)
maxl = - total(0.5*alog(2.*!dpi*err^2)) - chi/2.
return, maxl
end