newreg2.cpp
Go to the documentation of this file.
```00001 /*
00002  * \$Id: newreg2.cpp 1708 2014-02-28 21:43:41Z johnoel \$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <fvar.hpp>
00012
00017 dvariable robust_regression(const dvector& obs, const dvar_vector& pred,
00018   const dvariable& a)
00019 {
00020   if (obs.indexmin() != pred.indexmin() || obs.indexmax() != pred.indexmax() )
00021   {
00022     cerr << "Index limits on observed vector are not equal to the Index\n"
00023       "limits on the predicted vector in robust_reg_likelihood function\n";
00024   }
00025   RETURN_ARRAYS_INCREMENT(); //Need this statement because the function
00026                              //returns a variable type
00027   dvariable v_hat;
00028   double width=3.0;
00029   double pcon=0.05;
00030   double width2=width*width;
00031   dvariable a2;
00032   a2=a*a;
00033   dvar_vector diff = obs-pred;     // These are the residuals
00034   dvar_vector diff2 = pow(diff,2); // These are the squared residuals
00035   v_hat = mean(diff2)+1.e-80; // add 1.e-80 so that a perfect fit wont't
00036                               // produce log(0).
00037   double b=2.*pcon/(width*sqrt(PI));  // This is the weight for the
00038                                       // "robustifying" term
00039   dvariable log_likelihood = -sum(log((1.-pcon)*exp(-diff2/(2.*a2*v_hat))
00040     + b/(1.+pow(diff2/(width2*a2*v_hat),2))));
00041   log_likelihood += 0.5*diff.size()*log(a2*v_hat);
00042   RETURN_ARRAYS_DECREMENT(); // Need this to decrement the stack increment
00043                              // caused by RETURN_ARRAYS_INCREMENT();
00044   return(log_likelihood);
00045 }
00046
00051 dvariable robust_regression(const dvector& obs, const dvar_vector& pred,
00052   CGNU_DOUBLE a)
00053 {
00054   if (obs.indexmin() != pred.indexmin() || obs.indexmax() != pred.indexmax() )
00055   {
00056     cerr << "Index limits on observed vector are not equal to the Index\n"
00057       "limits on the predicted vector in robust_reg_likelihood function\n";
00058   }
00059   RETURN_ARRAYS_INCREMENT(); //Need this statement because the function
00060                              //returns a variable type
00061   dvariable v_hat;
00062   double width=3.0;
00063   double pcon=0.05;
00064   //double width2=width*width;
00065   dvariable a2;
00066   a2=a*a;
00067   dvar_vector diff = obs-pred;     // These are the residuals
00068   dvar_vector diff2 = square(diff); // These are the squared residuals
00069   v_hat = mean(diff2)+1.e-80; // add 1.e-80 so that a perfect fit wont't
00070                               // produce log(0).
00071   double b=2.*pcon/(width*sqrt(PI));  // This is the weight for the
00072                                       // "robustifying" term
00073   dvariable log_likelihood = -sum(log((1.-pcon)*exp(-diff2/(2.*a2*v_hat))
00074     + b/(1.+pow(diff2/(a2*v_hat),2))));
00075   log_likelihood += 0.5*diff.size()*log(a2*v_hat);
00076   RETURN_ARRAYS_DECREMENT(); // Need this to decrement the stack increment
00077                              // caused by RETURN_ARRAYS_INCREMENT();
00078   return(log_likelihood);
00079 }
```