mestim.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "libtpcmodel.h"
/*****************************************************************************/
 
/*****************************************************************************/
double mEstim(
  double *data,
  int nr,
  int iterNr,
  double cutoff
) {
  double theta, sum1, sum2, help;
 
  theta=dmedian(data, nr);
  for(int ii=0; ii<iterNr; ii++){
    sum1=sum2=0;
    for(int in=0; in<nr; in++){
      if(data[in]<0.9999*theta || data[in]>1.0001*theta){
        help=huber(data[in]-theta, cutoff);
        sum1=sum1+data[in]*help/(data[in]-theta);
        sum2=sum2+help/(data[in]-theta);
      } else {
        sum1=sum1+cutoff*data[in];
        sum2=sum2+cutoff; 
      }
    }
    theta=sum1/sum2;
  }
  return theta;
}
/*****************************************************************************/
 
/*****************************************************************************/
double huber(
  double x, 
  double b
) {
  double help;
 
  if(x<-b) help=-b; else help=x;
  if(help<b) {return help;} else {return b;}
}
/*****************************************************************************/
 
/*****************************************************************************/