Sympy

docs/requirements.txt

@@ -8,3 +8,4 @@ matplotlib
 nbinteract
 scipy
 jupyter-book
+sympy

ngboost/.vscode/settings.json

@@ -0,0 +1,3 @@
+{
+    "python.pythonPath": "D:\\src\\ngboost\\venv\\python.exe"
+}

ngboost/distns/__init__.py

@@ -9,3 +9,7 @@
 from .normal import Normal, NormalFixedVar  # NOQA
 from .poisson import Poisson  # NOQA
 from .t import T, TFixedDf, TFixedDfFixedVar  # NOQA
+from .gamma import Gamma
+from .beta import Beta
+from .logitnormal import LogitNormal
+from .betabernoulli import BetaBernoulli

ngboost/distns/beta.py

@@ -0,0 +1,83 @@
+import numpy as np
+import numpy
+import math
+import scipy as sp
+from scipy.stats import beta as dist
+from scipy.special import polygamma, beta
+
+from ngboost.distns.distn import RegressionDistn
+from ngboost.scores import CRPScore, LogScore
+
+import sympy as sym
+from sympy import stats as symstats
+from sympy.printing.lambdarepr import NumPyPrinter
+from sympy.utilities.lambdify import lambdastr
+
+logalpha, logbeta, x = sym.symbols('logalpha logbeta x')
+
+distr = symstats.Beta('dist', sym.exp(logalpha), sym.exp(logbeta))
+score = -sym.log(symstats.density( distr ).pdf(x))
+def neg_loglikelihood_sympy(logalpha, logbeta, x):
+    return score
+
+# Need to create str, because non-numpy functions are used: polygamma; and exp is not translated into numpy inside polygamma
+def newlambdify(args, funs):
+    funcstr = lambdastr(args, funs, printer=NumPyPrinter)
+    funcstr = funcstr.replace(
+        ' exp', 'numpy.exp'
+    )
+    return eval(funcstr)
+
+neg_loglikelihood = np.vectorize( newlambdify((logalpha, logbeta, x), neg_loglikelihood_sympy(logalpha, logbeta, x)) )
+D_0 = np.vectorize( newlambdify( (logalpha, logbeta, x), sym.diff(neg_loglikelihood_sympy(logalpha, logbeta, x), logalpha)) )
+D_1 = np.vectorize( newlambdify( (logalpha, logbeta, x), sym.diff(neg_loglikelihood_sympy(logalpha, logbeta, x), logbeta)) )
+# FI_0_0 = np.vectorize( newlambdify( (logalpha, logbeta), sym.factor(sym.expand(symstats.E(sym.factor(sym.expand(sym.diff(sym.diff(score, logalpha), logalpha))).subs(x, distr))))) )
+# FI_0_1 = np.vectorize( newlambdify( (logalpha, logbeta), sym.factor(sym.expand(symstats.E(sym.factor(sym.expand(sym.diff(sym.diff(score, logalpha), logbeta))).subs(x, distr))))) )
+# FI_1_0 = np.vectorize( newlambdify( (logalpha, logbeta), sym.factor(sym.expand(symstats.E(sym.factor(sym.expand(sym.diff(sym.diff(score, logbeta), logalpha))).subs(x, distr))))) )
+# FI_1_1 = np.vectorize( newlambdify( (logalpha, logbeta), sym.factor(sym.expand(symstats.E(sym.factor(sym.expand(sym.diff(sym.diff(score, logbeta), logbeta))).subs(x, distr))))) )
+
+
+
+class BetaLogScore(LogScore):
+
+    def score(self, Y):
+        return -dist.logpdf(Y, a=self.alpha, b=self.beta, loc=0, scale=1)
+        # return neg_loglikelihood(k=self.a, logtheta=self.logscale, x=Y)
+
+    def d_score(self, Y):
+        D = np.zeros((len(Y), 2)) # first col is dS/da, second col is dS/d(log(scale))
+        D[:, 0] = D_0(logalpha=self.logalpha, logbeta=self.logbeta, x=Y)
+        D[:, 1] = D_1(logalpha=self.logalpha, logbeta=self.logbeta, x=Y)
+        return D
+
+class Beta(RegressionDistn):
+
+    n_params = 2
+    scores = [BetaLogScore]
+
+    def __init__(self, params):
+        # save the parameters
+        self._params = params
+
+        # create other objects that will be useful later
+        self.logalpha = params[0]
+        self.logbeta = params[1]
+        self.alpha = np.exp(params[0])
+        self.beta = np.exp(params[1]) # since params[1] is log(scale)
+        self.dist = dist(a=self.alpha, b=self.beta)
+
+    def fit(Y):
+        alpha, beta, loc1, scale1  = dist.fit(Y, floc=0, fscale=1) # use scipy's implementation
+        return np.array([np.log(alpha), np.log(beta)])
+
+    def sample(self, m):
+        return np.array([self.dist.rvs() for i in range(m)])
+
+    def __getattr__(self, name): # gives us access to Laplace.mean() required for RegressionDist.predict()
+        if name in dir(self.dist):
+            return getattr(self.dist, name)
+        return None
+
+    @property
+    def params(self):
+        return {'alpha':self.alpha, 'beta':self.beta}

ngboost/distns/betabernoulli.py

@@ -0,0 +1,180 @@
+from scipy.stats import betabinom as dist
+import numpy as np
+from ngboost.distns.distn import RegressionDistn
+from ngboost.scores import LogScore
+from scipy.special import digamma, polygamma
+from array import array
+import sys
+
+class BetaBernoulliLogScore(LogScore):
+    def score(self, Y):
+        return -self.dist.logpmf(Y)
+
+    def d_alpha(self, Y):
+        return self.alpha * (
+            digamma(self.alpha + self.beta)
+            + digamma(Y + self.alpha)
+            - digamma(self.alpha + self.beta + 1)
+            - digamma(self.alpha)
+        )
+
+    def d_beta(self, Y):
+        return self.beta * (
+            digamma(self.alpha + self.beta)
+            + digamma(-Y + self.beta + 1)
+            - digamma(self.alpha + self.beta + 1)
+            - digamma(self.beta)
+        )
+
+    def d_alpha_alpha(self, Y):
+        return (
+            (polygamma(1, Y + self.alpha) +
+             polygamma(1, self.alpha + self.beta) -
+             polygamma(1, self.alpha + self.beta + 1) -
+             polygamma(1, self.alpha)
+             ) * self.alpha**2 + self.d_alpha(Y)
+        ) 
+
+    def d_alpha_beta(self, Y):
+        return (
+            self.alpha * self.beta * (
+                polygamma(1, self.alpha + self.beta) -
+                polygamma(1, self.alpha + self.beta + 1)
+            )
+        )
+
+    def d_beta_beta(self, Y):
+        return (
+            (polygamma(1, -Y + self.beta + 1) + 
+             polygamma(1, self.alpha + self.beta) - 
+             polygamma(1, self.alpha + self.beta + 1) - 
+             polygamma(1, self.beta)
+            ) * self.beta**2 + self.d_beta(Y)
+        )
+    def d_score(self, Y):
+        D = np.zeros(
+            (len(Y), 2)
+        )  # first col is dS/d(log(α)), second col is dS/d(log(β))
+        D[:, 0] = -self.d_alpha(Y)
+        D[:, 1] = -self.d_beta(Y)
+        return D
+
+    # Variance
+    def metric(self):
+        FI = np.zeros((self.alpha.shape[0], 2, 2))
+        FI[:, 0, 0] = (
+            self.d_alpha(0)*self.d_alpha(0)*self.dist.pmf(0) +
+            self.d_alpha(1)*self.d_alpha(1)*self.dist.pmf(1)
+        ) 
+        # FI[:, 1, 0] = (
+        #     self.d_alpha(0)*self.d_beta(0)*self.dist.pmf(0) +
+        #     self.d_alpha(1)*self.d_beta(1)*self.dist.pmf(1)
+        # ) 
+        # FI[:, 0, 1] = (
+        #     self.d_alpha(0)*self.d_beta(0)*self.dist.pmf(0) +
+        #     self.d_alpha(1)*self.d_beta(1)*self.dist.pmf(1)
+        # ) 
+        FI[:, 1, 1] = (
+            self.d_beta(0)*self.d_beta(0)*self.dist.pmf(0) +
+            self.d_beta(1)*self.d_beta(1)*self.dist.pmf(1)
+        ) 
+        return FI
+
+    # Hessian
+    # def metric(self):
+    #     FI = np.zeros((self.alpha.shape[0], 2, 2))
+    #     FI[:, 0, 0] = self.d_alpha_alpha(0) * self.dist.pmf(0) + self.d_alpha_alpha(1) * self.dist.pmf(1)
+    #     # FI[:, 1, 0] = self.d_alpha_beta(0) * self.dist.pmf(0) + self.d_alpha_beta(1) * self.dist.pmf(1)
+    #     # FI[:, 0, 1] = self.d_alpha_beta(0) * self.dist.pmf(0) + self.d_alpha_beta(1) * self.dist.pmf(1)
+    #     FI[:, 1, 1] = self.d_beta_beta(0) * self.dist.pmf(0) + self.d_beta_beta(1) * self.dist.pmf(1)
+    #     return FI
+
+class BetaBernoulli(RegressionDistn):
+
+    n_params = 2
+    scores = [BetaBernoulliLogScore]
+
+    def __init__(self, params):
+        # save the parameters
+        self._params = params
+
+        # create other objects that will be useful later
+        self.log_alpha = params[0]
+        self.log_beta = params[1]
+        self.alpha = np.exp(self.log_alpha)
+        self.beta = np.exp(self.log_beta)
+        self.dist = dist(n=1, a=self.alpha, b=self.beta)
+
+    def fit(Y):
+        def fit_alpha_beta_py(success, alpha0=1.5, beta0=5, niter=1000):
+            # based on https://github.com/lfiaschi/fastbetabino/blob/master/fastbetabino.pyx
+
+            # This optimisation works for Beta-Binomial distribution in general. 
+            # For Beta-Bernoulli it's simplified by fixing the trials to 1.
+            trials = np.ones_like(Y)
+
+            alpha_old = alpha0
+            beta_old = beta0
+
+            for it in range(niter):
+
+                alpha = (
+                    alpha_old
+                    * (
+                        sum(
+                            digamma(c + alpha_old) - digamma(alpha_old)
+                            for c, i in zip(success, trials)
+                        )
+                    )
+                    / (
+                        sum(
+                            digamma(i + alpha_old + beta_old)
+                            - digamma(alpha_old + beta_old)
+                            for c, i in zip(success, trials)
+                        )
+                    )
+                )
+
+                beta = (
+                    beta_old
+                    * (
+                        sum(
+                            digamma(i - c + beta_old) - digamma(beta_old)
+                            for c, i in zip(success, trials)
+                        )
+                    )
+                    / (
+                        sum(
+                            digamma(i + alpha_old + beta_old)
+                            - digamma(alpha_old + beta_old)
+                            for c, i in zip(success, trials)
+                        )
+                    )
+                )
+
+                # print('alpha {} | {}  beta {} | {}'.format(alpha,alpha_old,beta,beta_old))
+                sys.stdout.flush()
+
+                if np.abs(alpha - alpha_old) and np.abs(beta - beta_old) < 1e-10:
+                    # print('early stop')
+                    break
+
+                alpha_old = alpha
+                beta_old = beta
+
+            return alpha, beta
+
+        alpha, beta = fit_alpha_beta_py(Y)
+        return np.array([np.log(alpha), np.log(beta)])
+
+    def sample(self, m):
+        return np.array([self.dist.rvs() for i in range(m)])
+
+    def __getattr__(self, name):  
+        if name in dir(self.dist):
+            return getattr(self.dist, name)
+        return None
+
+    @property
+    def params(self):
+        return {"alpha": self.alpha, "beta": self.beta}