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?