Add initial proportional priors

conjugate/distributions.py

@@ -42,6 +42,7 @@
 import numpy as np
 
 from scipy import stats, __version__ as scipy_version
+from scipy.special import gammaln
 
 from conjugate._typing import NUMERIC
 from conjugate.plot import (
@@ -671,3 +672,119 @@ class CompoundGamma(ContinuousPlotDistMixin, SliceMixin):
     alpha: NUMERIC
     beta: NUMERIC
     lam: NUMERIC
+
+
+@dataclass
+class GammaKnownRateProportional:
+    """Gamma known rate proportional distribution.
+
+    Args:
+        a: prod of observations
+        b: number of observations
+        c: number of observations
+
+    """
+
+    a: NUMERIC
+    b: NUMERIC
+    c: NUMERIC
+
+    def approx_log_likelihood(
+        self, alpha: NUMERIC, beta: NUMERIC, ln=np.log, gammaln=gammaln
+    ):
+        """Approximate log likelihood.
+
+        Args:
+            alpha: shape parameter
+            beta: known rate parameter
+            ln: log function
+            gammaln: log gamma function
+
+        Returns:
+            log likelihood up to a constant
+
+        """
+        return (
+            (alpha - 1) * ln(self.a)
+            + alpha * self.c * ln(beta)
+            - self.b * gammaln(alpha)
+        )
+
+
+@dataclass
+class GammaProportional:
+    """Gamma proportional distribution.
+
+    Args:
+        p: product of r observations
+        q: sum of s observations
+        r: number of observations for p
+        s: number of observations for q
+
+    """
+
+    p: NUMERIC
+    q: NUMERIC
+    r: NUMERIC
+    s: NUMERIC
+
+    def approx_log_likelihood(
+        self, alpha: NUMERIC, beta: NUMERIC, ln=np.log, gammaln=gammaln
+    ):
+        """Approximate log likelihood.
+
+        Args:
+            alpha: shape parameter
+            beta: rate parameter
+            ln: log function
+            gammaln: log gamma function
+
+        Returns:
+            log likelihood up to a constant
+
+        """
+        return (
+            (alpha - 1) * ln(self.p)
+            - self.q * beta
+            - self.r * gammaln(alpha)
+            + self.s * alpha * ln(beta)
+        )
+
+
+@dataclass
+class BetaProportional:
+    """Beta proportional distribution.
+
+    Args:
+        p: product of observations
+        q: product of complements
+        k: number of observations
+
+    """
+
+    p: NUMERIC
+    q: NUMERIC
+    k: NUMERIC
+
+    def approx_log_likelihood(
+        self, alpha: NUMERIC, beta: NUMERIC, ln=np.log, gammaln=gammaln
+    ):
+        """Approximate log likelihood.
+
+        Args:
+            alpha: shape parameter
+            beta: shape parameter
+            ln: log function
+            gammaln: log gamma function
+
+        Returns:
+            log likelihood up to a constant
+
+        """
+        return (
+            self.k * gammaln(alpha + beta)
+            + alpha * ln(self.p)
+            + beta * ln(self.q)
+            - self.k * gammaln(alpha)
+            - self.k * gammaln(beta)
+        )

conjugate/models.py

@@ -9,10 +9,13 @@
 
 from conjugate.distributions import (
     Beta,
+    BetaProportional,
     CompoundGamma,
     Dirichlet,
     DirichletMultinomial,
     Gamma,
+    GammaKnownRateProportional,
+    GammaProportional,
     NegativeBinomial,
     BetaNegativeBinomial,
     BetaBinomial,
@@ -1164,3 +1167,84 @@ def pareto_gamma(
     beta_post = gamma_prior.beta + ln_x_total - n * ln(x_m)
 
     return Gamma(alpha=alpha_post, beta=beta_post)
+
+
+def gamma(
+    x_total: NUMERIC,
+    x_prod: NUMERIC,
+    n: NUMERIC,
+    proportional_prior: GammaProportional,
+) -> GammaProportional:
+    """Posterior distribution for a gamma likelihood.
+
+    Inference on alpha and beta
+
+    Args:
+        x_total: sum of all outcomes
+        x_prod: product of all outcomes
+        n: total number of samples in x_total and x_prod
+        proportional_prior: GammaProportional prior
+
+    Returns:
+        GammaProportional posterior distribution
+
+    """
+    p_post = proportional_prior.p * x_prod
+    q_post = proportional_prior.q + x_total
+    r_post = proportional_prior.r + n
+    s_post = proportional_prior.s + n
+
+    return GammaProportional(p=p_post, q=q_post, r=r_post, s=s_post)
+
+
+def gamma_known_rate(
+    x_prod: NUMERIC,
+    n: NUMERIC,
+    beta: NUMERIC,
+    proportional_prior: GammaKnownRateProportional,
+) -> GammaKnownRateProportional:
+    """Posterior distribution for a gamma likelihood.
+
+    The rate beta is assumed to be known.
+
+    Args:
+        x_prod: product of all outcomes
+        n: total number of samples in x_prod
+        beta: known rate parameter
+
+    Returns:
+        GammaKnownRateProportional posterior distribution
+
+    """
+    a_post = proportional_prior.a * x_prod
+    b_post = proportional_prior.b + n
+    c_post = proportional_prior.c + n
+
+    return GammaKnownRateProportional(a=a_post, b=b_post, c=c_post)
+
+
+def beta(
+    x_prod: NUMERIC,
+    one_minus_x_prod: NUMERIC,
+    n: NUMERIC,
+    proportional_prior: BetaProportional,
+) -> BetaProportional:
+    """Posterior distribution for a Beta likelihood.
+
+    Inference on alpha and beta
+
+    Args:
+        x_prod: product of all outcomes
+        one_minus_x_prod: product of all (1 - outcomes)
+        n: total number of samples in x_prod and one_minus_x_prod
+        proportional_prior: BetaProportional prior
+
+    Returns:
+        BetaProportional posterior distribution
+
+    """
+    p_post = proportional_prior.p * x_prod
+    q_post = proportional_prior.q * one_minus_x_prod
+    k_post = proportional_prior.k + n
+
+    return BetaProportional(p=p_post, q=q_post, k=k_post)

tests/test_models.py

@@ -10,6 +10,7 @@
 
 from conjugate.distributions import (
     Beta,
+    BetaProportional,
     CompoundGamma,
     Dirichlet,
     Pareto,
@@ -21,6 +22,8 @@
     Normal,
     StudentT,
     MultivariateStudentT,
+    GammaProportional,
+    GammaKnownRateProportional,
 )
 from conjugate.models import (
     get_binomial_beta_posterior_params,
@@ -43,6 +46,9 @@
     normal_normal_inverse_gamma_posterior_predictive,
     gamma_known_shape,
     gamma_known_shape_posterior_predictive,
+    gamma,
+    gamma_known_rate,
+    beta,
 )
 
 rng = np.random.default_rng(42)
@@ -418,3 +424,59 @@ def test_gamma_known_shape(shape) -> None:
         alpha=shape, gamma=posterior
     )
     assert isinstance(posterior_predictive, CompoundGamma)
+
+
+def test_gamma_proportional_model() -> None:
+    true_alpha, true_beta = 2, 6
+    true = Gamma(alpha=true_alpha, beta=true_beta)
+
+    n_samples = 15
+    samples = true.dist.rvs(size=n_samples, random_state=0)
+
+    prior = GammaProportional(1, 1, 1, 1)
+    posterior = gamma(
+        x_total=samples.sum(),
+        x_prod=np.prod(samples),
+        n=n_samples,
+        proportional_prior=prior,
+    )
+
+    assert isinstance(posterior, GammaProportional)
+
+    # TODO: add a comparison to the true for prior and posterior
+
+
+def test_gamma_known_rate() -> None:
+    true_alpha, true_beta = 2, 6
+    true = Gamma(alpha=true_alpha, beta=true_beta)
+
+    n_samples = 15
+    samples = true.dist.rvs(size=n_samples, random_state=0)
+
+    prior = GammaKnownRateProportional(1, 1, 1)
+    posterior = gamma_known_rate(
+        x_prod=np.prod(samples),
+        n=n_samples,
+        beta=true_beta,
+        proportional_prior=prior,
+    )
+
+    assert isinstance(posterior, GammaKnownRateProportional)
+
+
+def test_beta_proportional_model() -> None:
+    true_alpha, true_beta = 2, 6
+    true = Beta(alpha=true_alpha, beta=true_beta)
+
+    n_samples = 15
+    samples = true.dist.rvs(size=n_samples, random_state=0)
+
+    prior = BetaProportional(0.25, 0.25, 1)
+    posterior = beta(
+        x_prod=np.prod(samples),
+        one_minus_x_prod=np.prod(1 - samples),
+        n=n_samples,
+        proportional_prior=prior,
+    )
+
+    assert isinstance(posterior, BetaProportional)