Dear all,
Can someone fit two function to a data and perform Akaike information criterion test on them? A sample code would be great!

For pedagogical purposes, I would post the code that I have written yesterday based on the code that was written by davekw7x. For AIC formula readers are referred to:
http://www.nr.com/forum/newreply.php?do=newreply&noquote=1&p=5309

//
// Driver for akaike
//
//Amir Vahid + davekw7x

//
// Always include nr3.h first
//
#include ".../code/nr3.h"

// Will use a Gaussian random number generator
// to simulate measurement noise. Include
// these to be able to declare a Normaldev object
#include ".../code/ran.h"
#include ".../code/gamma.h"
#include ".../code//deviates.h"

// The following are needed for Fitsvd
#include ".../code/svd.h"
#include ".../code/fitsvd.h"

//
// A global variable used in the fpoly() function.
//
// Note that you can't simply #include unmodified fit_examples.h
// for this problem since the distribution file sets fpoly_np
// equal to 10, and we will be using a fourth degree polynomial
//
// fpoly_np is set to the number of coefficients of a fourth-degree
// polynomial.
//
Int fpoly_np = 5;
//
// Here's the prototype of fpoly:
//
VecDoub fpoly(const Doub x);


int main(void)
{
Int i;

Int numpoints = 201;
Int npol = fpoly_np;

VecDoub x(numpoints); // Independent variable values
VecDoub y(numpoints); // "Measured" values
VecDoub sig(numpoints); // Estimates of standard deviations of measured values

// A "random" object to generate Gaussian noise is initialized with
// a seed obtained from current time
//
Doub perturb = 0.02; // 2% perturbation of "measured" value
Normaldev normaldev(0.0, perturb, time(0));
//
// For test purposes, might want to use a fixed seed
//Normaldev normaldev(0.0, perturb, 876543321);
//
//

cout << fixed << setprecision(6);
//
// Create a set of data points for a polynomial
//
// The polynomial is 1 + 2x + 3x^2 + 4x^3 + 5x^4
//
// Add some "noise" to the values
//
Doub xstart = 0.0;
Doub xmax = 2.0;
Doub deltax = (xmax - xstart) / (numpoints - 1.0);
Doub xx;
Doub dev;
Doub AIC;
Doub num_Parm;

num_Parm = npol + 1;
for (i = 0, xx = xstart; i < numpoints; i++, xx += deltax) {
x[i] = xx;
y[i] = 1.0 + xx * (2.0 + xx * (3.0 + xx * (4.0 + xx * 5.0)));
dev = normaldev.dev();
y[i] += y[i] * dev;

}
//
// For perturb = 0.02, we estimate "measurement error" for all points
// to be 2% of the measured value.
//
for (i = 0; i < numpoints; i++) {
sig[i] = perturb * y[i];
}

//
// Invoke the constructor to bind the data and the function
//
Fitsvd fitsvd(x, y, sig, fpoly);


fitsvd.fit();
double rss=0;
//cout << "poly coeff are:\n";
for (int j=0; j<npol; j ++) cout <<fitsvd.a[j]<< "\t";
//cout << endl;
for (int i=0; i<numpoints; i++){
double y_fit=0;
for (int j=0; j<npol; j ++){
//cout <<x[i]<< endl;
//cout << pow( x[i], j)<< endl;
y_fit += fitsvd.a[j] * pow( x[i], j) ;

}
//cout << y_fit << endl;
rss += SQR( y[i] - y_fit);
}


cout << endl << "Polynomial coefficients:" << endl << endl;
cout << setw(11) << "Actual"
<< setw(13) << "Estimated"
<< setw(17) << "Err Estimate"
<< endl;

for (i = 0; i < npol; i++) {
cout << setw(12) << (Doub)(i + 1);
cout << setw(12) << fitsvd.a[i] << " +-";
cout << setw(11) << sqrt(fitsvd.covar[i][i]) << endl;
}
if ((numpoints/num_Parm) < 40){
AIC = numpoints * log(rss/numpoints) + 2 * num_Parm + 2 * num_Parm * (num_Parm + 1) / (numpoints - num_Parm - 1);
}
if ((numpoints/num_Parm) >= 40){
AIC = numpoints * log(rss/numpoints) + 2 * num_Parm;
}

cout << endl << "Number of points = " << numpoints
<< ", Chi-squared = " << fitsvd.chisq << ", RSS values = " << rss <<", AIC values = "<< AIC << endl;
return 0;
}

VecDoub fpoly(const Doub x) {
Int j;
VecDoub p(fpoly_np);
p[0]=1.0;
for (j=1;j<fpoly_np;j++) p[j]=p[j-1]*x;
return p;
}