aic.c File Reference

Routines for model selection and weighting using Akaike's information criteria. More...

#include "libtpcmodel.h"

Go to the source code of this file.

Functions
double	aicSS (double ss, const int n, const int k)
int	parFreeNr (const int n, double *pLower, double *pUpper)
	Calculate the number of free parameters.
int	aicWeights (double *aic, double *w, int n)
double	aicWeightedAvg (double *w, double *p, int n)
double	aicModel (double *w, int n)

Detailed Description

Routines for model selection and weighting using Akaike's information criteria.

https://www.turkupetcentre.net/reports/tpcmod0016.pdf

Author

Kaisa Sederholm, Vesa Oikonen

Definition in file aic.c.

Function Documentation

aicModel()

double aicModel	(	double *	w,
		int	n )

Calculates a value describing the relative goodness of models, based on an array of model weights.

See also

aicSS, aicWeightedAvg, aicWeights

Returns

Returns the weighted average of model number.

Parameters

w	Array of weights
n	Length of array

Definition at line 132 of file aic.c.

  {
  int i;
  double avg;
 
  if(n<1 || w==NULL) avg=0.0;
  else for(i=0, avg=0.0; i<n; i++) avg+=w[i]*(double)(i+1);
  return(avg);
}

aicSS()

double aicSS	(	double	ss,
		const int	n,
		const int	k )

Computation of AICc in the special case of sum-of-squares optimization from the SS, nr of fitted points and nr of fitted parameters.

If variance is different between the data points, weighted SS must be given.

See also

parFreeNr, dftWSampleNr, aicWeights, aicWeightedAvg

Returns

Returns the AIC value.

Parameters

ss	Sum-of-Squares of the fit
n	Sample size, i.e. nr of fitted data points
k	Number of fitted model parameters; AICc calculation is valid only when (n-k)>1.

Definition at line 20 of file aic.c.

  {
  if(!(ss>=0.0) || n<1 || k<0 || (n-k)<2) return(nan(""));
 
  double aic=0.0, bias_adj, css;
  int dr, dv;
  
  dr=n-k-1; dv=2*k*(k+1);
  if(dr>0) bias_adj=((double)dv)/((double)dr); else bias_adj=0.0;
  if(ss<1.0e-50) css=1.0e-50; else css=ss; /* Because log(0) is an error */
  if(n>0) aic= n*log(css/(double)n) + 2.0*(double)k + bias_adj;
  return(aic);
}

Referenced by imgNoiseTemplate().

aicWeightedAvg()

double aicWeightedAvg	(	double *	w,
		double *	p,
		int	n )

Computation of the Akaike weighted model parameter average. Requires arrays of AIC weight values, and corresponding parameter values.

See also

aicSS, aicModel, aicWeights

Returns

Returns the weighted average.

Parameters

w	Array of weights
p	Array of parameters
n	Lengths of arrays

Definition at line 108 of file aic.c.

  {
  int i;
  double avg;
 
  /* Check data */
  if(n<1 || w==NULL || p==NULL) return(1);
  for(i=0, avg=0.0; i<n; i++) avg+=w[i]*p[i];
  return(avg);
}

aicWeights()

int aicWeights	(	double *	aic,
		double *	w,
		int	n )

Computation of the Akaike weights for model averaging. Requires an array of AIC values, and an output array for weights.

See also

aicSS, aicWeightedAvg, aicModel

Returns

Returns 0, if OK.

Parameters

aic	Array of AICs
w	Array of weights (output)
n	Lengths of arrays

Definition at line 74 of file aic.c.

  {
  int i, mini;
  double minaic, sume;
 
  /* Check data */
  if(n<1 || aic==NULL || w==NULL) return(1);
  if(n==1) {w[0]=1.0; return(0);}
  /* Find out which model gave the smallest AIC */
  mini=0; for(i=1; i<n; i++) if(aic[i]<aic[mini]) mini=i;
  minaic=aic[mini];
  /* Compute relative weights for each model */
  for(i=0, sume=0.0; i<n; i++) {
    w[i]=exp(-0.5*(aic[i]-minaic));
    sume+=w[i];
  }
  if(sume==0.0) return(2);
  for(i=0; i<n; i++) w[i]/=sume;
  return(0);
}

parFreeNr()

int parFreeNr	(	const int	n,
		double *	pLower,
		double *	pUpper )

Calculate the number of free parameters.

Model parameters can be fixed by setting lower and upper limit to equal values. This function simply checks the limits for each parameter.

Returns

Returns the number of free parameters.

See also

aicSS, dftWSampleNr

Parameters

n	Nr of parameters
pLower	Lower limits (array of length n)
pUpper	Upper limits (array of length n)

Definition at line 50 of file aic.c.

  {
  if(n<1 || pLower==NULL || pUpper==NULL) return(0);
  int nf=0;
  for(int i=0; i<n; i++) {
    double range=pUpper[i]-pLower[i];
    if(range>1.0E-10) nf++;
  }
  return(nf);
}