Python scipy.stats.norm.logpdf() Examples

The following are 30 code examples of scipy.stats.norm.logpdf(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module scipy.stats.norm , or try the search function .
Example #1
Source File: linear.py    From cgpm with Apache License 2.0 6 votes vote down vote up
def logpdf_marginal(self, z):
        return norm.logpdf(z, scale=1) 
Example #2
Source File: linear.py    From cgpm with Apache License 2.0 6 votes vote down vote up
def logpdf_joint(self, x, y):
        return multivariate_normal.logpdf(
            np.array([x,y]), np.array([0,0]),
            np.array([[1,1-self.noise],[1-self.noise,1]])) 
Example #3
Source File: test_comparisons.py    From ChainConsumer with MIT License 6 votes vote down vote up
def test_aic_posterior_dependence():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    p2 = norm.logpdf(d, scale=2)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=1000)
    c.add_chain(d, posterior=p2, num_free_params=1, num_eff_data_points=1000)
    aics = c.comparison.aic()
    assert len(aics) == 2
    assert aics[0] == 0
    expected = 2 * np.log(2)
    assert np.isclose(aics[1], expected, atol=1e-3) 
Example #4
Source File: test_likelihoods.py    From revrand with Apache License 2.0 6 votes vote down vote up
def test_gaussian():

    # Test we can at match a Gaussian distribution from scipy

    mu = 0
    var = 2
    dist = lk.Gaussian()

    x = np.linspace(-10, 10, 100)

    p1 = norm.logpdf(x, loc=mu, scale=np.sqrt(var))
    p2 = dist.loglike(x, mu, var)

    np.allclose(p1, p2)

    p1 = norm.cdf(x, loc=mu, scale=np.sqrt(var))
    p2 = dist.cdf(x, mu, var)

    np.allclose(p1, p2) 
Example #5
Source File: mvknn.py    From cgpm with Apache License 2.0 6 votes vote down vote up
def logpdf(self, rowid, targets, constraints=None, inputs=None):
        constraints = self.populate_constraints(rowid, targets, constraints)
        # XXX Disable logpdf queries without constraints.
        if inputs:
            raise ValueError('Prohibited inputs: %s' % (inputs,))
        if not constraints:
            raise ValueError('Provide at least one constraint: %s'
                % (constraints,))
        self._validate_simulate_logpdf(rowid, targets, constraints)
        # Retrieve the dataset and neighborhoods.
        dataset, neighborhoods = self._find_neighborhoods(targets, constraints)
        models = [self._create_local_model_joint(targets, dataset[n])
            for n in neighborhoods]
        # Compute logpdf in each neighborhood and simple average.
        lp = [m.logpdf(targets) for m in models]
        return gu.logsumexp(lp) - np.log(len(models)) 
Example #6
Source File: run_regression.py    From Doubly-Stochastic-DGP with Apache License 2.0 6 votes vote down vote up
def _event_handler(self, manager):
        minibatch_size = 1000
        S = 100
        means, vars = [], []
        for mb in range(-(-len(Xs) // minibatch_size)):
            m, v = model.predict_y(Xs[mb * minibatch_size:(mb + 1) * minibatch_size, :], S)
            means.append(m)
            vars.append(v)
        mean_SND = np.concatenate(means, 1)
        var_SDN = np.concatenate(vars, 1)

        mean_ND = np.average(mean_SND, 0)

        test_err = np.average(Y_std * np.mean((Ys - mean_ND) ** 2.0) ** 0.5)
        test_nll_ND = logsumexp(norm.logpdf(Ys * Y_std, mean_SND * Y_std, var_SDN ** 0.5 * Y_std), 0, b=1 / float(S))
        test_nll = np.average(test_nll_ND)

        summary, step = model.enquire_session().run([self.summary, global_step],
                                                    feed_dict={self._full_test_err: test_err,
                                                               self._full_test_nlpp: test_nll})
        self.file_writer.add_summary(summary, step) 
Example #7
Source File: test_analysis.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_max_likelihood_4(self):
        c = ChainConsumer()
        data = np.random.multivariate_normal([0, 0], [[1, 0], [0, 1]], size=10000)
        posterior = norm.logpdf(data).sum(axis=1)
        data[:, 1] += 2
        c.add_chain(data, parameters=["x", "y"], posterior=posterior, name="A")
        c.add_chain(data + 3, parameters=["x", "y"], name="B")
        result = c.analysis.get_max_posteriors(parameters="x", chains="A", squeeze=False)
        assert len(result) == 1
        x = result[0]["x"]
        assert np.isclose(x, 0, atol=0.05) 
Example #8
Source File: test_comparisons.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_bic_0():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=1000)
    bics = c.comparison.bic()
    assert len(bics) == 1
    assert bics[0] == 0 
Example #9
Source File: test_comparisons.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_bic_posterior_dependence():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    p2 = norm.logpdf(d, scale=2)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=1000)
    c.add_chain(d, posterior=p2, num_free_params=1, num_eff_data_points=1000)
    bics = c.comparison.bic()
    assert len(bics) == 2
    assert bics[0] == 0
    expected = 2 * np.log(2)
    assert np.isclose(bics[1], expected, atol=1e-3) 
Example #10
Source File: test_comparisons.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_bic_parameter_dependence():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=1000)
    c.add_chain(d, posterior=p, num_free_params=2, num_eff_data_points=1000)
    bics = c.comparison.bic()
    assert len(bics) == 2
    assert bics[0] == 0
    expected = np.log(1000)
    assert np.isclose(bics[1], expected, atol=1e-3) 
Example #11
Source File: test_comparisons.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_bic_data_dependence2():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=2, num_eff_data_points=1000)
    c.add_chain(d, posterior=p, num_free_params=3, num_eff_data_points=500)
    bics = c.comparison.bic()
    assert len(bics) == 2
    assert bics[0] == 0
    expected = 3 * np.log(500) - 2 * np.log(1000)
    assert np.isclose(bics[1], expected, atol=1e-3) 
Example #12
Source File: test_comparisons.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_dic_0():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p)
    dics = c.comparison.dic()
    assert len(dics) == 1
    assert dics[0] == 0 
Example #13
Source File: test_comparisons.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_dic_posterior_dependence():
    d = norm.rvs(size=1000000)
    p = norm.logpdf(d)
    p2 = norm.logpdf(d, scale=2)
    c = ChainConsumer()
    c.add_chain(d, posterior=p)
    c.add_chain(d, posterior=p2)
    bics = c.comparison.dic()
    assert len(bics) == 2
    assert bics[1] == 0
    dic1 = 2 * np.mean(-2 * p) + 2 * norm.logpdf(0)
    dic2 = 2 * np.mean(-2 * p2) + 2 * norm.logpdf(0, scale=2)
    assert np.isclose(bics[0], dic1 - dic2, atol=1e-3) 
Example #14
Source File: test_analysis.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_max_likelihood_1(self):
        c = ChainConsumer()
        data = np.random.multivariate_normal([0, 0], [[1, 0], [0, 1]], size=10000)
        posterior = norm.logpdf(data).sum(axis=1)
        data[:, 1] += 1
        c.add_chain(data, parameters=["x", "y"], posterior=posterior, name="A")
        result = c.analysis.get_max_posteriors()
        x, y = result["x"], result["y"]
        assert np.isclose(x, 0, atol=0.05)
        assert np.isclose(y, 1, atol=0.05) 
Example #15
Source File: test_analysis.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_max_likelihood_2(self):
        c = ChainConsumer()
        data = np.random.multivariate_normal([0, 0], [[1, 0], [0, 1]], size=10000)
        posterior = norm.logpdf(data).sum(axis=1)
        data[:, 1] += 2
        c.add_chain(data, parameters=["x", "y"], posterior=posterior, name="A")
        c.add_chain(data + 3, parameters=["x", "y"], name="B")
        result = c.analysis.get_max_posteriors(parameters=["x", "y"], chains="A")
        x, y = result["x"], result["y"]
        assert np.isclose(x, 0, atol=0.05)
        assert np.isclose(y, 2, atol=0.05) 
Example #16
Source File: multivariate_normal.py    From ngboost with Apache License 2.0 5 votes vote down vote up
def nll(self, Y):
        try:
            E = Y["Event"]
            T = np.log(Y["Time"] + eps)

            (
                mu0_given_1,
                mu1_given_0,
                var0_given_1,
                var1_given_0,
            ) = self.conditional_dist(T)

            cens = (1 - E) * (
                dist.logpdf(T, loc=self.loc[:, 1], scale=self.cov[:, 1, 1] ** 0.5)
                + np.log(
                    eps + 1 - dist.cdf(T, loc=mu0_given_1, scale=var0_given_1 ** 0.5)
                )
            )
            uncens = E * (
                dist.logpdf(T, loc=self.loc[:, 0], scale=self.cov[:, 0, 0] ** 0.5)
                + np.log(
                    eps + 1 - dist.cdf(T, loc=mu1_given_0, scale=var1_given_0 ** 0.5)
                )
            )
            return -(cens + uncens)
        except:
            diff = Y - self.loc
            M = diff[:, None, :] @ self.cov_inv @ diff[:, :, None]
            half_log_det = np.log(eps + np.diagonal(self.L, axis1=1, axis2=2)).sum(-1)
            const = self.p / 2 * np.log(eps + 2 * np.pi)
            logpdf = -const - half_log_det - 0.5 * M.flatten()
            return -logpdf 
Example #17
Source File: test_sample.py    From pyPESTO with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def gaussian_llh(x):
    return float(norm.logpdf(x)) 
Example #18
Source File: test_plotter.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_plotter_extents6(self):
        c = ChainConsumer()
        for mid in np.linspace(-1, 1, 3):
            data = np.random.normal(loc=0, size=1000)
            posterior = norm.logpdf(data)
            data += mid
            c.add_chain(data, parameters=['x'], posterior=posterior, plot_point=True, plot_contour=False)

        c.configure()
        minv, maxv = c.plotter._get_parameter_extents("x", c.chains)
        assert np.isclose(minv, -1, atol=0.01)
        assert np.isclose(maxv, 1, atol=0.01) 
Example #19
Source File: train.py    From focus with GNU General Public License v3.0 5 votes vote down vote up
def _fit_cis_herit(y, K_cis, X=None, compute_lrt=True):
    log = logging.getLogger(pyfocus.LOG)

    try:
        from glimix_core.lmm import LMM
        from numpy_sugar.linalg import economic_qs_linear
    except ImportError as ie:
        log.error("Training submodule requires glimix-core>=2.0.0 and numpy-sugar to be installed.")
        raise

    from scipy.stats import norm
    from scipy.linalg import lstsq

    if X is None:
        X = np.ones((len(y), 1))

    K_cis = economic_qs_linear(K_cis)
    lmm = LMM(y, X, K_cis)
    lmm.fit(verbose=False)

    fixed_betas = lmm.beta
    logl_1 = lmm.lml()

    cis_scale = lmm.v0
    noise_scale = lmm.v1
    fe_scale = lmm.fixed_effects_variance

    if compute_lrt:
        n, p = X.shape
        # reduced model is just OLS regression for fixed-effects
        fixed_betas_0, sosqs, ranks, svals = lstsq(X, y)
        s2e = sosqs / len(y)  # LMM also uses MLE estimation, so don't adjust for bias right now

        logl_0 = np.sum(norm.logpdf(y, loc=np.dot(X, fixed_betas_0), scale=np.sqrt(s2e)))
        pval = _lrt_pvalue(logl_0, logl_1)
        log.debug("Estimated cis-h2g = {} (P = {})".format(cis_scale / (cis_scale + noise_scale + fe_scale), pval))
    else:
        pval = None
        log.debug("Estimated cis-h2g = {}".format(cis_scale / (cis_scale + noise_scale + fe_scale)))

    return fe_scale, cis_scale, noise_scale, logl_1, fixed_betas, pval 
Example #20
Source File: normal_test.py    From MXFusion with Apache License 2.0 5 votes vote down vote up
def test_log_pdf(self, dtype, mean, mean_isSamples, var, var_isSamples,
                     rv, rv_isSamples, num_samples):
        from scipy.stats import norm

        isSamples_any = any([mean_isSamples, var_isSamples, rv_isSamples])
        rv_shape = rv.shape[1:] if rv_isSamples else rv.shape
        n_dim = 1 + len(rv.shape) if isSamples_any and not rv_isSamples else len(rv.shape)
        mean_np = numpy_array_reshape(mean, mean_isSamples, n_dim)
        var_np = numpy_array_reshape(var, var_isSamples, n_dim)
        rv_np = numpy_array_reshape(rv, rv_isSamples, n_dim)
        log_pdf_np = norm.logpdf(rv_np, mean_np, np.sqrt(var_np))
        normal = Normal.define_variable(shape=rv_shape, dtype=dtype).factor
        mean_mx = mx.nd.array(mean, dtype=dtype)
        if not mean_isSamples:
            mean_mx = add_sample_dimension(mx.nd, mean_mx)
        var_mx = mx.nd.array(var, dtype=dtype)
        if not var_isSamples:
            var_mx = add_sample_dimension(mx.nd, var_mx)
        rv_mx = mx.nd.array(rv, dtype=dtype)
        if not rv_isSamples:
            rv_mx = add_sample_dimension(mx.nd, rv_mx)
        variables = {normal.mean.uuid: mean_mx, normal.variance.uuid: var_mx, normal.random_variable.uuid: rv_mx}
        log_pdf_rt = normal.log_pdf(F=mx.nd, variables=variables)

        assert np.issubdtype(log_pdf_rt.dtype, dtype)
        assert array_has_samples(mx.nd, log_pdf_rt) == isSamples_any
        if isSamples_any:
            assert get_num_samples(mx.nd, log_pdf_rt) == num_samples
        if np.issubdtype(dtype, np.float64):
            rtol, atol = 1e-7, 1e-10
        else:
            rtol, atol = 1e-4, 1e-5
        assert np.allclose(log_pdf_np, log_pdf_rt.asnumpy(), rtol=rtol, atol=atol) 
Example #21
Source File: normal_mean_precision_test.py    From MXFusion with Apache License 2.0 5 votes vote down vote up
def test_log_pdf(self, dtype, mean, mean_is_samples, precision, precision_is_samples,
                     rv, rv_is_samples, num_samples):
        is_samples_any = any([mean_is_samples, precision_is_samples, rv_is_samples])
        rv_shape = rv.shape[1:] if rv_is_samples else rv.shape
        n_dim = 1 + len(rv.shape) if is_samples_any and not rv_is_samples else len(rv.shape)
        mean_np = numpy_array_reshape(mean, mean_is_samples, n_dim)
        precision_np = numpy_array_reshape(precision, precision_is_samples, n_dim)
        rv_np = numpy_array_reshape(rv, rv_is_samples, n_dim)
        log_pdf_np = norm.logpdf(rv_np, mean_np, np.power(precision_np, -0.5))

        var = NormalMeanPrecision.define_variable(shape=rv_shape, dtype=dtype).factor
        mean_mx = mx.nd.array(mean, dtype=dtype)
        if not mean_is_samples:
            mean_mx = add_sample_dimension(mx.nd, mean_mx)
        precision_mx = mx.nd.array(precision, dtype=dtype)
        if not precision_is_samples:
            precision_mx = add_sample_dimension(mx.nd, precision_mx)
        rv_mx = mx.nd.array(rv, dtype=dtype)
        if not rv_is_samples:
            rv_mx = add_sample_dimension(mx.nd, rv_mx)
        variables = {var.mean.uuid: mean_mx, var.precision.uuid: precision_mx, var.random_variable.uuid: rv_mx}
        log_pdf_rt = var.log_pdf(F=mx.nd, variables=variables)

        assert np.issubdtype(log_pdf_rt.dtype, dtype)
        assert array_has_samples(mx.nd, log_pdf_rt) == is_samples_any
        if is_samples_any:
            assert get_num_samples(mx.nd, log_pdf_rt) == num_samples
        if np.issubdtype(dtype, np.float64):
            rtol, atol = 1e-7, 1e-10
        else:
            rtol, atol = 1e-4, 1e-5
        assert np.allclose(log_pdf_np, log_pdf_rt.asnumpy(), rtol=rtol, atol=atol) 
Example #22
Source File: numpy_backend.py    From pyhf with Apache License 2.0 5 votes vote down vote up
def normal_logpdf(self, x, mu, sigma):
        # this is much faster than
        # norm.logpdf(x, loc=mu, scale=sigma)
        # https://codereview.stackexchange.com/questions/69718/fastest-computation-of-n-likelihoods-on-normal-distributions
        root2 = np.sqrt(2)
        root2pi = np.sqrt(2 * np.pi)
        prefactor = -np.log(sigma * root2pi)
        summand = -np.square(np.divide((x - mu), (root2 * sigma)))
        return prefactor + summand

    # def normal_logpdf(self, x, mu, sigma):
    #     return norm.logpdf(x, loc=mu, scale=sigma) 
Example #23
Source File: svd.py    From graspy with Apache License 2.0 5 votes vote down vote up
def _compute_likelihood(arr):
    """
    Computes the log likelihoods based on normal distribution given
    a 1d-array of sorted values. If the input has no variance,
    the likelihood will be nan.
    """
    n_elements = len(arr)
    likelihoods = np.zeros(n_elements)

    for idx in range(1, n_elements + 1):
        # split into two samples
        s1 = arr[:idx]
        s2 = arr[idx:]

        # deal with when input only has 2 elements
        if (s1.size == 1) & (s2.size == 1):
            likelihoods[idx - 1] = -np.inf
            continue

        # compute means
        mu1 = np.mean(s1)
        if s2.size != 0:
            mu2 = np.mean(s2)
        else:
            # Prevent numpy warning for taking mean of empty array
            mu2 = -np.inf

        # compute pooled variance
        variance = ((np.sum((s1 - mu1) ** 2) + np.sum((s2 - mu2) ** 2))) / (
            n_elements - 1 - (idx < n_elements)
        )
        std = np.sqrt(variance)

        # compute log likelihoods
        likelihoods[idx - 1] = np.sum(norm.logpdf(s1, loc=mu1, scale=std)) + np.sum(
            norm.logpdf(s2, loc=mu2, scale=std)
        )

    return likelihoods 
Example #24
Source File: shared.py    From widedeepnetworks with Apache License 2.0 5 votes vote down vote up
def baseline_log_density(Y_test, train_mean=0., train_std=1.):
    #Evaluate the log density associated with 
    #a base line model that uses a single normal density with
    #the mean and std of the training data.
    log_pdfs = norm.logpdf(Y_test, loc=train_mean, scale=train_std)
    return log_pdfs.mean() 
Example #25
Source File: test_comparisons.py    From ChainConsumer with MIT License 5 votes vote down vote up
def test_bic_fail_no_data_points():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1)
    bics = c.comparison.bic()
    assert len(bics) == 1
    assert bics[0] is None 
Example #26
Source File: linear.py    From cgpm with Apache License 2.0 5 votes vote down vote up
def logpdf_conditional(self, w, z):
        mean = self.conditional_mean(z)
        var = self.conditional_variance(z)
        return norm.logpdf(w, loc=mean, scale=np.sqrt(var)) 
Example #27
Source File: dots.py    From cgpm with Apache License 2.0 5 votes vote down vote up
def logpdf_joint(self, x, y):
        return gu.logsumexp([np.log(.25)
                + norm.logpdf(x, loc=mx, scale=self.noise)
                + norm.logpdf(y, loc=my, scale=self.noise)
            for (mx,my) in zip(self.mx, self.my)]) 
Example #28
Source File: dots.py    From cgpm with Apache License 2.0 5 votes vote down vote up
def logpdf_marginal(self, z):
        return gu.logsumexp(
            [np.log(.5) + norm.logpdf(z, loc=mx, scale=self.noise)
            for mx in set(self.mx)]) 
Example #29
Source File: normal_trunc.py    From cgpm with Apache License 2.0 5 votes vote down vote up
def calc_predictive_logp(x, mu, sigma, l, h):
        return norm.logpdf(x, loc=mu, scale=sigma) 
Example #30
Source File: normal_trunc.py    From cgpm with Apache License 2.0 5 votes vote down vote up
def calc_log_prior(mu, sigma, alpha, beta, l, h):
        log_sigma = gamma.logpdf(sigma, alpha, scale=beta)
        log_mu = uniform.logpdf(mu, loc=l, scale=l+h)
        return log_mu + log_sigma