aic.c

Go to the documentation of this file.

 
/****************************************************************************/
#include "libtpcmodel.h"
/****************************************************************************/
 
/****************************************************************************/
double aicSS(
  double ss,
  const int n,
  const int k
) {
  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);
}
/****************************************************************************/
 
/****************************************************************************/
int parFreeNr(
  const int n,
  double *pLower,
  double *pUpper
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int aicWeights(
  double *aic,
  double *w,
  int n
) {
  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);
}
/****************************************************************************/
 
/****************************************************************************/
double aicWeightedAvg(
  double *w,
  double *p,
  int n
) {
  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);
}
/****************************************************************************/
 
/****************************************************************************/
double aicModel(
  double *w,
  int n
) {
  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);
}
/****************************************************************************/
 
/****************************************************************************/