Does this implementation make sense?
class multivariate_skewnorm(multi_rv_generic):
def __init__(self, mu, shape, cov=None):
self.dim = len(shape)
self.shape = np.asarray(shape)
self.mean = np.asarray(mu)
self.cov = np.eye(self.dim) if cov is None else np.asarray(cov)
def pdf(self, x):
return np.exp(self.logpdf(x))
def logpdf(self, x):
x = mvn._process_quantiles(value, self.dim)
pdf = mvn(self.mean, self.cov).logpdf(x)
cdf = norm(0, 1).logcdf(np.dot(x, self.shape))
return _squeeze_output(np.log(2) + pdf + cdf)
def rvs(self, size=1, random_state=43):
random_state = self._get_random_state(random_state)
aCa = self.shape @ self.cov @ self.shape
delta = (1 / np.sqrt(1 + aCa)) * self.cov @ self.shape
cov_star = np.block([[np.ones(1), delta],
[delta[:, None], self.cov]])
if size == None:
size = 1
else:
size = size[0]
x = random_state.multivariate_normal(np.zeros(self.dim+1), cov_star, size=size)
x0, x1 = x[:, 0], x[:, 1:]
inds = x0 <= 0
x1[inds] = -1 * x1[inds]
return x1
@as_op(itypes=[pt.dmatrix, pt.dmatrix, pt.dmatrix, pt.dvector], otypes=[pt.dscalar])
def logp(value, mu, cov, alpha):
return multivariate_skewnorm(mu, alpha, cov).logpdf(value)
def random(mu, cov, alpha, rng=None, size=None):
return multivariate_skewnorm(mu, alpha, cov).rvs(size=size, random_state=23)