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
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
Cheers,
Régis
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
