code: def gaussfunc(y, ybar, sigma):
      
      x = old_div((y - ybar), (np.sqrt(2.) * sigma))
      t = old_div(1.0, (1.0 + .3275911 * abs(x)))
      erf = 1.0 - np.exp(-x * x) * t * (.254829592 - t * (.284496736 -
                                                          t * (1.421413741 - t * (1.453152027 - t * 1.061405429))))
      erf = abs(erf)
      sign = old_div(x, abs(x))
      return 0.5 * (1.0 + sign * erf)
docs: cumulative normal distribution function of the variable y
      with mean ybar,standard deviation sigma
      uses expression 7.1.26 from Abramowitz & Stegun
      accuracy better than 1.5e-7 absolute
      Parameters
      _________
      y : input variable
      ybar : mean
      sigma : standard deviation