Kent distribution

Python usage
1

Hi all at OT,

I need to implement the Kent Distribution Kent distribution - Wikipedia, more precisely a method to sample from it. Is there a “standard” way to achieve this in OT, meaning using existent Python class that would make things easier in a “full OT world” ?

Best regards

Ti

2

Hi Ti,

The difficulty is that this distribution does not have a density function wrt the Lebesgue measure in R^3, but wrt the uniform distribution over the unit sphere S^2. It makes it difficult to sample using generic techniques. The trick is to use a 2D parameterization of the unit sphere, then to sample a companion distribution in the parameter space, and translate the parameter sample into a 3D sample.

Something like this should do the trick (not sure it is 100% correct but you will see the spirit)

```python
import openturns as ot
import openturns.experimental as otexp
import math

class KentDistribution(ot.PythonDistribution):
    def __init__(self, kappa=1.0, beta=0.25, gamma1=[1.0, 0.0, 0.0], gamma2=[0.0, 1.0, 0.0], gamma3=[0.0, 0.0, 1.0], epsilon=1e-15):
        super().__init__(3)
        # Here you add the tests on the arguments
        # 0 <= 2*beta < kappa
        # gamma1, gamma2, gamma3 3D unit vectors
        # pairwise orthogonals
        self.kappa_ = kappa
        self.beta_ = beta
        self.gamma1_ = ot.Point(gamma1)
        self.gamma2_ = ot.Point(gamma2)
        self.gamma3_ = ot.Point(gamma3)
        self.epsilon_ = epsilon
        # Compute the normalization factor
        factor = 0.0
        go = True
        j = 0
        logBeta = math.log(beta)
        logHalfKappa = math.log(0.5 * kappa)
        while go:
            delta = math.exp(ot.SpecFunc.LogGamma(j + 0.5) - ot.SpecFunc.LogGamma(j + 1.0) + 2.0 * j * logBeta - (2.0 * j + 0.5) * logHalfKappa + ot.SpecFunc.LogBesselInu(kappa, 2.0 * j + 0.5))
            factor += delta
            go = delta < epsilon * factor
        self.normalization_ = factor
        def logUnscaledPDFPy(X):
            return [math.log(self.computePDFParam(X))]
        logUnscaledPDF = ot.PythonFunction(2, 1, logUnscaledPDFPy)
        isScaled = True
        self.sampler_ = otexp.RatioOfUniforms(logUnscaledPDF, ot.Interval([0.0]*2, [math.pi, 2.0 * math.pi]), isScaled)
        
    def getRange(self):
        return ot.Interval([-1.0]*3, [1.0]*3)

    def getRealization(self):
        theta, phi = self.sampler_.getRealization()
        sTheta = math.sin(theta)
        cTheta = math.cos(theta)
        sPhi = math.sin(phi)
        cPhi = math.cos(phi)
        return [sTheta * cPhi, sTheta * sPhi, cTheta]

    def computeCDF(self, X):
        raise

    def computePDFParam(self, X):
        theta, phi = X
        sTheta = math.sin(theta)
        cTheta = math.cos(theta)
        sPhi = math.sin(phi)
        cPhi = math.cos(phi)
        return sTheta * self.computePDF([sTheta * cPhi, sTheta * sPhi, cTheta])

    def computePDF(self, X):
        # Assume that X is a dimension 3 vector
        X = ot.Point(X)
        normX = X.norm()
        if abs(normX - 1.0) > self.epsilon_:
            return 0.0
        return math.exp(self.kappa_ * X.dot(self.gamma1_) + self.beta_ * (X.dot(self.gamma2_)**2 - X.dot(self.gamma3_)**2)) / self.normalization_

if __name__ == '__main__':
    # Helper to build orthonormal bases
    def buildGammas():
        M = ot.SquareMatrix(3, ot.DistFunc.rNormal(9))
        Q, R = M.computeQR()
        e1 = ot.Point([1.0, 0.0, 0.0])
        e2 = ot.Point([0.0, 1.0, 0.0])
        e3 = ot.Point([0.0, 0.0, 1.0])
        return [Q*e1, Q*e2, Q*e3]

    # kappa, beta, gamma1, gamma2, gamma3
    params = list()
    params.append([2.5, 0.5, [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])
    params.append([5.0, 2.0] + buildGammas())
    params.append([10.0, 3.0] + buildGammas())
    print(params)
    import matplotlib.pyplot as plt

    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')

    for p in params:
        print(p)
        kappa, beta, gamma1, gamma2, gamma3 = p
        distribution = ot.Distribution(KentDistribution(kappa, beta, gamma1, gamma2, gamma3))
        sample = distribution.getSample(500)
        ax.scatter(sample[:,0], sample[:, 1], sample[:, 2], marker='o')

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')

    plt.show()


It produces

Figure_1
Figure_1640×480 83.5 KB

Cheers,

Régis

3

Hi Regis, thank you very much ! as simple as this !!! what would be the next steps (if any) for full use of the Kent Distribution in ot for UQ ou SA tasks ? can I simply combine it with other distributions ???

thansk again, Best Regards, Thierry

4

Hi Thierry,
You can find the documentation of PythonDistribution here: https://openturns.github.io/openturns/latest/user_manual/\_generated/openturns.PythonDistribution.html
and an example on how to do it here:

You will find additional information on how to create and manipulate your own distributions.
A nice thing is that you asked your question while I was implementing the vonMises-Fisher distribution (see Added the Von Mises Fisher distribution and factory by regislebrun · Pull Request #3220 · openturns/openturns · GitHub). As the Kent distribution is closely related to this one, I implemented it too (see Added Kent distribution by regislebrun · Pull Request #3223 · openturns/openturns · GitHub). It is still a work in progress, but it is almost completed and you should have it in the experimental module in the next release (autumn). You will be able to use it sooner, thanks to the nightly builds. Wait for a week after the PR is merged.

Cheers

Régis