x,y,sigma = var("x,y,sigma")
xn,yn=var("x_n,y_n")
Yn=var("Y_n")
G,a0,gamma,R=var("G_n a_0 gamma R")
yb=function("yb",nargs=2,latex_name="\overline{Y_n}")
def L(yb):
return -log(sqrt(2*pi)*sigma)-(1/2)*((Yn-yb(x,y))/sigma)**2
def j(expr,sym):
return [expr.diff(i) for i in sym]
J= j(L(yb),[x,y,sigma])