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))))))