Implement transformed_dist

yli/distributions.py

@@ -243,7 +243,43 @@ class betaoddsrat_gen(stats.rv_continuous):
 # See beta_oddsratio.cdf
 BETA_ODDSRATIO_CDF_TERMS = [100]
 
-beta_oddsratio = betaoddsrat_gen(name='beta_oddsratio', a=0, shapes='a1,b1,a2,b2')  # a = lower bound of support
+beta_oddsratio = betaoddsrat_gen(name='beta_oddsratio', a=0)  # a = lower bound of support
+
+class transformed_gen(stats.rv_continuous):
+	"""
+	Represents a transformation, Y, of a "base" random variable, X, where f(Y) = X
+	Hence Y.pdf(x) = X.pdf(f(x)) * f'(x)
+	
+	e.g. if X is a model parameter, then transformed_dist(X, f=np.exp, fprime=np.exp, finv=np.log) is distribution of the log-parameter
+	"""
+	
+	def _pdf(self, x, base, f, fprime, finv):
+		return base.pdf(np.vectorize(f)(x)) * np.vectorize(fprime)(x)
+	
+	def _cdf(self, x, base, f, fprime, finv):
+		return base.cdf(np.vectorize(f)(x))
+	
+	def _ppf(self, p, base, f, fprime, finv):
+		return np.vectorize(finv)(base.ppf(p))
+	
+	# Call implementations directly to bypass SciPy args checking, which chokes on "base" parameter
+	
+	def pdf(self, x, base, f, fprime, finv):
+		if np.isscalar(x):
+			return np.atleast_1d(self._pdf(x, base, f, fprime, finv))[0]
+		return self._pdf(x, base, f, fprime, finv)
+	
+	def cdf(self, x, base, f, fprime, finv):
+		if np.isscalar(x):
+			return np.atleast_1d(self._cdf(x, base, f, fprime, finv))[0]
+		return self._cdf(x, base, f, fprime, finv)
+	
+	def ppf(self, x, base, f, fprime, finv):
+		if np.isscalar(x):
+			return np.atleast_1d(self._ppf(x, base, f, fprime, finv))[0]
+		return self._ppf(x, base, f, fprime, finv)
+
+transformed_dist = transformed_gen(name='transformed')
 
 ConfidenceInterval = collections.namedtuple('ConfidenceInterval', ['lower', 'upper'])