Exercise 28 solutions

# start from ipython --pylab

from scipy.optimize import minimize

# our data
x = array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype="float")
y = array([18, 5, 17, 38, 40, 106, 188, 234, 344, 484], dtype="float")

# plot the data
figure()
plot(x,y,'b.',markersize=15)
xlabel("X")
ylabel("Y")

# define our cost function
def e28costfun(b,x,y,ff):
        y_est = ff(b,x)
        return sum((y_est-y)**2)

# define the model of the data
def e28mod(b,x):
        return b[0] + b[1]*x + b[2]*(x**2) + b[3]*(x**3)

# initial guess at beta parameters
b0 = rand(4,1)
# optimize!
opt_result = minimize(e28costfun, b0, method="Nelder-Mead", args=(x,y,e28mod))
print opt_result

# get predicted values of y
b = opt_result["x"]
c = opt_result["fun"]
xp = linspace(min(x),max(x),100)
yp = e28mod(b, xp)

# plot them
plot(xp, yp, 'r-', linewidth=2)

x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
y = [18, 5, 17, 38, 40, 106, 188, 234, 344, 484];

% the "@" symbol denotes a function handle
b = fminsearch(@e28cost, randn(1,4), [], x, y, @e28mod);

plot(x, y, 'b.', 'markersize', 30)
hold on

xhat = linspace(min(x),max(x),200);
yhat = e28mod(b,xhat);
plot(xhat, yhat, 'g-', 'linewidth', 2)

function out = e28cost(b,x,y,ff)

        y_est = ff(b,x);
        out = sum((y_est-y).^2);

function out = e28mod(b,x)

        out = b(1) + b(2).*x + b(3).*(x.^2) + b(4).*(x.^3);

# define the cost function
e28cost <- function(b,x,y,fmod) {
	yhat <- fmod(b,x)
	J <- sum((yhat-y)**2)
	J
}

# define the data model
e28mod <- function(b,x) {
	yhat <- b[1] + b[2]*x + b[3]*x*x + b[4]*x*x*x
	yhat
}

# the data
x <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y <- c(18, 5, 17, 38, 40, 106, 188, 234, 344, 484)

# plot x vs y
plot(x, y, "p", col="blue", pch=19, cex=1.5, xlab="X", ylab="Y")

# initial guess
b0 <- runif(4)

# optimize!
optim_res <- optim(par=b0, fn=e28cost, method="Nelder-Mead", gr=NULL, x, y, e28mod)
print(optim_res)
b_opt <- optim_res$par
min_cost <- optim_res$value
cat(paste("initial guess was", b0, "\n"))
cat(paste("the optimal value of b is", b_opt, "\n"))
cat(paste("the minimum cost found is", min_cost, "\n"))

# plot line of best fit
xhat <- seq(min(x), max(x), (max(x)-min(x))/100)
yhat <- b_opt[1] + b_opt[2]*xhat + b_opt[3]*xhat*xhat + b_opt[4]*xhat*xhat*xhat
lines(xhat, yhat, col="red", lwd=2)

// e28.c
// compile with: gcc -o e28 e28.c nmsimplex.c

#include <stdio.h>
#include "nmsimplex.h"

// define a struct that will hold extra stuff to be passed
// to the cost function
typedef struct {
  int n;            // length of x,y arrays
  double *x, *y;    // the x,y data arrays
                    // below: a pointer to a function that computes
                    // y estimates based on b, x and n
  double *(*modfn)(double b[], double x[], int n);
} optim_extras;

// our cost function
double e28cost(double b[], void *extras) {
  optim_extras *stuff = (optim_extras *)extras;
  double J = 0.0;
  double *x = stuff->x;
  double *y = stuff->y;
  int i;
  double *(*fn)() = stuff->modfn;                 // the model function
  double *yhat = fn(b, x, stuff->n);              // compute yhat
  for (i=0; i<stuff->n; i++) {
    J = J + ( (yhat[i]-y[i]) * (yhat[i]-y[i]) );  // compute cost
  }
  free(yhat);
  return J;
}

// our model function to compute yhat from b, x and n
double *e28mod(double b[], double x[], int n) {
  double *yhat = malloc(n * sizeof(double));
  int i;
  for (i=0; i<n; i++) {
    yhat[i] = b[0] + (b[1]*x[i]) + (b[2]*x[i]*x[i]) + (b[3]*x[i]*x[i]*x[i]);
  }
  return yhat;
}

// the main program
//
int main(int argc, char *argv[]) {

  double x[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
  double y[] = {18, 5, 17, 38, 40, 106, 188, 234, 344, 484};

  optim_extras stuff;     // declare struct of extra stuff
  stuff.n = 10;           // x,y array length
  stuff.x = x;            // x data array
  stuff.y = y;            // y data array
  stuff.modfn = &e28mod;  // the address of the model function e28mod()

  double b[] = {.1,.1,.1,.1};                                   // our initial guess
  double f_min = simplex(e28cost, b, 4, 1e-8, 1, NULL, &stuff); // optimize!
  printf("f_min=%f\n", f_min);
  printf("x_min=(%f,%f,%f,%f)\n", b[0],b[1],b[2],b[3]);

  return 0;
}