import scipy.optimize as op
#
# Model
#
def rvmodel(x,ampl,offset,phase):
    return offset+ampl*np.cos(2*np.pi*x+phase)

#
# Likelihood function
#
def lnlike(theta, x, y, yerr):
    ampl, offset, phase = theta
    model = rvmodel(x,ampl,offset,phase)
    sigma2 = yerr**2
    return -0.5*(np.sum((y-model)**2*/sigma2 + np.log(2*np.pi*sigma2)))

#
# Bounds
#
bounds = [(10.,100.),(-50.,50.),(0.,2*np.pi)]

#
# Maximum likelihood fit
#
nll = lambda *args: -lnlike(*args)
result = op.minimize(nll, [50., 10., 3.], args=(x, y, yerr), bounds=bounds)
ampl_ml, offset_ml, phase_ml = result["x"]