code: def _func(a, p, r, v):
      
      # eq. 2.3
      f = a[0]*math.log(r-1.) + \
          a[1]*math.log(r-1.)**2 + \
          a[2]*math.log(r-1.)**3 + \
          a[3]*math.log(r-1.)**4
  
      # eq. 2.7 and 2.8 corrections
      if r == 3:
          f += -0.002 / (1. + 12. * _phi(p)**2)
  
          if v <= 4.364:
              f += 1./517. - 1./(312.*(v,1e38)[np.isinf(v)])
          else:
              f += 1./(191.*(v,1e38)[np.isinf(v)])
  
      return -f
docs: calculates f-hat for the coefficients in a, probability p,
      sample mean difference r, and degrees of freedom v.