diff --git a/yli/distributions.py b/yli/distributions.py index 5c5e518..0c8312a 100644 --- a/yli/distributions.py +++ b/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'])