Cauchy distribution in OT?

Python usage
1

Hi!
I would like to perform a Monte-Carlo simulation which fails to converge. To create it, I selected the Cauchy distribution https://en.wikipedia.org/wiki/Cauchy_distribution.
However, I could not find it in the library: did I miss something?
Best regards,
Michaël
PS
As a workaround, I created it based on the PythonDistribution class:

class Cauchy(ot.PythonDistribution):
    def __init__(self, location=0.0, scale=1.0):
        super(Cauchy, self).__init__(1)
        self.location = location
        self.scale = scale

    def computeCDF(self, x):
        u = x[0]
        xi = (u - self.location) / self.scale
        p = np.arctan(xi) / np.pi + 0.5
        return p

    def computePDF(self, x):
        u = x[0]
        xi = (u - self.location) / self.scale
        y = 1.0 / (np.pi * self.scale * (1.0 + xi ** 2))
        return y

    def getRange(self):
        epsilon = ot.ResourceMap.GetAsScalar("Distribution-DefaultCDFEpsilon")
        z_min = np.tan(np.pi * (epsilon - 0.5))
        x_min = self.location + self.scale * z_min
        z_max = np.tan(np.pi * (0.5 - epsilon))
        x_max = self.location + self.scale * z_max
        return ot.Interval(x_min, x_max)

    def computeQuantile(self, prob, tail=False):
        if tail:
            z = np.tan(np.pi * (0.5 - prob))
        else:
            z = np.tan(np.pi * (prob - 0.5))
        x = self.location + self.scale * z
        return [x]

However, the getSample method is quite slow … until the computeQuantile method is implemented!

2

Did you noticed that the Cauchy distribution is the Student distribution with nu=1?

4

If you look at StudentFunction.cxx, you see

  if (nu == 1.0)
  {
    Scalar value = 0.0;
    if (std::abs(u) < 0.025373628595705897178)
    {
      const Scalar u2 = u * u;
      value = (-0.31830988618379067153 + (1.0471975511965977462 + (0.68902837067332933726 + (0.64766070854027820799 + 0.63921549794217821540 * u2) * u2) * u2) * u2) / u;
    }
    else value = std::tan((u - 0.5) * M_PI);
    return (tail == (p < 0.5) ? -value : value);
  }

so you get an improved speed and accuracy for small values of prob