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

## No comments:

## Post a Comment