# Creating a Custom Hypersecant distribution

I have been working on creating a custom distribution for Hypersecant distribution. However, I am not sure if I am doing the custom class correctly. I need some advice because I have been try to create a Bayesian Fusion approach using the custom distribution, but, I always keep on getting different results when I re-run my code.

So, I appreciate some advice on the design of the custom Hypersecant distribution. Also, what could be the reason for getting different distribution when I re-run my code of Bayesian Fusion?

``````class Hypsecant(Continuous):
"""
Parameters
----------
mu: float
Location parameter.
sigma: float
Scale parameter (sigma > 0). Converges to the standard deviation as nu increases. (only required if lam is not specified)
lam: float
Scale parameter (lam > 0). Converges to the precision as nu increases. (only required if sigma is not specified)
"""

def __init__(self, mu=None, sigma=None, lam=None, sd=None, *args, **kwargs):
# super().__init__(*args, **kwargs)
# super(Hypsecant, self).__init__(*args, **kwargs)
if sd is not None:
sigma = sd
warnings.warn("sd is deprecated, use sigma instead", DeprecationWarning)
lam, sigma = get_tau_sigma(tau=lam, sigma=sigma)
self.lam = lam = tt.as_tensor_variable(lam)
self.sigma = self.sd = sigma = tt.as_tensor_variable(sigma)
self.mean = self.median = self.mode = self.mu = mu = tt.as_tensor_variable(mu)
self.variance = tt.as_tensor_variable(1)

assert_negative_support(mu, 'mu', 'Hypsecant')
assert_negative_support(sigma, 'sigma (lam)', 'Hypsecant')

# return super(Hypsecant, self).__init__(shape=[mu, sigma], *args, **kwargs)
return super().__init__(*args, **kwargs)

def random(self, point=None, size=None):
"""
Draw random values from Hypsecant distribution.
Parameters
----------
point: dict, optional
Dict of variable values on which random values are to be
conditioned (uses default point if not specified).
size: int, optional
Desired size of random sample (returns one sample if not
specified).
Returns
-------
array
"""
# mu = self.mu
# lam = self.lam
# sigma = self.sigma
# # sd = self.sd
# mu = self.mu
# lam, sigma = get_tau_sigma(sigma=sigma)

mu, sigma = draw_values([self.mu, self.sigma], point=point, size=size)
return generate_samples(stats.hypsecant.rvs, loc=mu, scale=sigma, dist_shape=self.shape, size=size)

def logp(self, value):
"""
Calculate log-probability of Hypsecant distribution at specified value.
Parameters
----------
value: numeric
Value(s) for which log-probability is calculated. If the log probabilities for multiple
values are desired the values must be provided in a numpy array or theano tensor
Returns
-------
TensorVariable
"""

# mu = self.mu
# lam = self.lam
sigma = self.sigma
# # sd = self.sd
mu = self.mu
# lam, sigma = get_tau_sigma(sigma=sigma)

# Px = pm.math.log(1.0/(math.pi * pm.math.cosh(value)))
# Px = pm.math.log(1.0/(math.pi * pm.math.cosh(((value - mu) / sigma))))
# Px = pm.math.log((1.0/math.pi) * (1 / pm.math.cosh(((value - mu) / sigma))))
Px = pm.math.log((1.0/math.pi) * (1 / pm.math.cosh((math.pi / 2) * value)))
# Px = pm.math.log((1.0/(2.0 * sigma)) * (1 / (pm.math.cosh((math.pi / 2) * ((value - mu) / sigma)))))
# Px = pm.math.log((1.0/(2.0 * sigma * ((pm.math.cosh((math.pi / 2) * ((value - mu) / sigma)))))))
return bound(Px)
# return bound(Px, mu, lam, sigma)

def logcdf(self, value):
"""
Compute the log of the cumulative distribution function for Hypsecant distribution
at the specified value.
Parameters
----------
value: numeric
Value(s) for which log CDF is calculated. If the log CDF for multiple
values are desired the values must be provided in a numpy array or theano tensor.
Returns
-------
TensorVariable
"""

# # mu = self.mu
# # lam = self.lam
# # # self.sigma = self.sd = sigma = tt.as_tensor_variable(sigma)
# mu = self.mu
# lam = self.lam
sigma = self.sigma
# # sd = self.sd
mu = self.mu
# lam, sigma = get_tau_sigma(sigma=sigma)

# return bound(pm.math.log(2.0 / math.pi * tt.math.atan(tt.math.exp(value))), mu, lam, sigma)
# return bound(pm.math.log(2.0 / math.pi * tt.math.atan(tt.math.exp(value))))
# return bound(pm.math.log(2.0 / math.pi * tt.math.atan(tt.math.exp((math.pi / 2.0) * ((value - mu) / sigma)))))
# return bound(pm.math.log((2.0 / math.pi) * (tt.math.atan(tt.math.exp((math.pi / 2.0) * ((value - mu) / sigma))))))
# return bound(pm.math.log((2.0 / math.pi) * (tt.math.atan(tt.math.exp((math.pi / 2.0) * ((value - mu) / sigma))))))
return bound(pm.math.log((2.0 / math.pi) * (tt.math.atan(tt.math.exp((math.pi / 2.0) * ((value)))))))
# return bound((2.0 / math.pi) * (tt.math.atan(tt.math.exp((math.pi / 2.0) * ((value - mu) / sigma)))))
# return bound(pm.math.log((2.0 / math.pi) * (tt.math.atan(tt.math.exp((math.pi / 2.0) * (value))))))
# return bound(pm.math.log((2.0 / math.pi * tt.math.atan(tt.math.exp(((value - mu) / sigma))))))
# return bound(pm.math.log(2.0 / math.pi * (tt.math.atan(tt.math.exp((math.pi / 2.0) * (value))))))
``````

Could it be because of the `def random`?

If you want to jump to PyMC V4 there’s a developer guide for new distributions: Implementing a Distribution — PyMC 0+untagged.310.g244c37d.dirty documentation

If you just need something simple you can also use DensityDist, which in V4 also allows for a random method (I don’t recall if it did in V3): pymc.DensityDist — PyMC 0+untagged.310.g244c37d.dirty documentation

Otherwise it is not clear what your problem is from your post alone.

Thank you very much for the reply.

To understand further about the custom distribution class configuration:

• What are the necessary components needed to design the class custom distribution?
• Is it necessary to include the `def random` in the custom distribution?
• What is the reason for getting different fused distributions when implementing the bayesian fusion part? Could it be because of `def random` code section in the custom distribution? For clarification about the fusion part, please find the following question: How to create Bayesian data fusion in python with pymc3?.