C:Optimization F:OPTIM L:function myfunc(w,z) L:{ L:x=w[1];y=w[2] L:t=x^2+(y-2)^2+exp(x+y)-x-y L:return(t) L:} L: L:function gradient(x,fun,Z) L:/* L:Calculate gradient of fun numerically L: L:Required arguments: L:x: Numerical values of variables L:fun: Character string, name of the optimized function without parentheses L:Z: List of arguments of fun L: L:Value: L:Numerical vector of gradient at x. L:*/ L: L:{ L:b=x L:k=nrows(x) L:eps=0.00001 L:zz=seq(0,0,count=k) L:for (i=1,k) L:{ L:b[i]=b[i]-eps L:z1=parse(fun+"(b,Z)") L:b[i]=b[i]+2*eps L:z2=parse(fun+"(b,Z)") L:b[i]=b[i]-eps L:zz[i]=(z2-z1)/(2*eps) L:} L:return(zz) L:} L: L:function hessian(x,fun,Z) L:/* L:Calculate gradient of fun numerically L: L:Required arguments: L:x: Numerical values of variables L:fun: Character string, name of the optimized function without parentheses L:Z: List of arguments of fun L: L:Value: L:Numerical matrix of Hessian at x. L:*/ L: L:{ L:b=x L:k=nrows(x) L:eps1=0.000001 L:eps2=0.00001 L:zz=unit(k) L:for (i=1,k) L:{ L:for(j=1,k) L:{ L:b[j]=b[j]-eps2 L:b[i]=b[i]-eps1 L:z1=parse(fun+"(b,Z)") L:b[i]=b[i]+2*eps2 L:z2=parse(fun+"(b,Z)") L:b[i]=b[i]-eps2 L:zz1=(z2-z1)/(2*eps2) L:b[j]=b[j]+2*eps1 L:b[i]=b[i]-eps2 L:z1=parse(fun+"(b,Z)") L:b[i]=b[i]+2*eps2 L:z2=parse(fun+"(b,Z)") L:b[i]=b[i]-eps2 L:zz2=(z2-z1)/(2*eps2) L:b[j]=b[j]-eps1 L:zz[i,j]=(zz2-zz1)/(2*eps1) L:} L:} L:return(zz) L:} L: L:function optimize(x,fun,Z) L:/* L:Unconstrained multivariate Gauss-Newton optimization L: L:Required arguments: L:x: Numerical vector of initial x-values L:fun: Character string, name of the optimized function without parentheses L: L:Value: L:A list with components: L:$x: Optimal values of x L:$obj: Value of FUN at extreme L:$iter: Number of iterations performed L:$grad: Value of gradient at optimum L:$hess: Hessian matrix of second derivatives at optimum L: L:Example: L:x=vec(1,1) L:Z=list(x=x[1],y=x[2]) L:fun="myfunc" L:a=optimize(x,"2*myfunc",Z) L:a$x // Optimal values of x L:a$obj // Value of FUN at extreme L:eigenval(a$hess) // Eigenvalues of hessian at extreme (all negative = maximum, all positive = minimum) L: L:*/ L:{ L:maxiter=100 L:x0=x L:obj0=parse(fun+"(x,Z)") L:gg=gradient(x,fun,Z) L:gnorm=norm(gg) L:iter=0 L:while(or( and(gt(gnorm,0.000001) , lt(iter,maxiter) ),not(iter)) ) L:{ L:hh=hessian(x,fun,Z) L:gg=gradient(x,fun,Z) L:step=-inv(hh)#gg L:x=x+step L:obj=parse(fun+"(x,Z)") L:gnorm=norm(gg) L:iter=iter+1 L:} // end While L:return(list(x=x, obj=obj, iter=iter,grad=gg,hess=hh)) L:} L: L: