Python scipy.stats.norm 模块，logpdf() 实例源码

def ln_prior(p, joker_params):
    P, phi0, ecc, omega, s, K, *v_terms = p

    lnp = 0.

    # TODO: more repeated code here and hard-coded priors
    if ecc < 0 or ecc > 1:
        return -np.inf

    lnp += beta.logpdf(ecc, 0.867, 3.03) # Kipping et al. 2013

    # TODO: do we need P_min, P_max here?

    if not joker_params._fixed_jitter:
        # DFM's idea: wide, Gaussian prior in log(s^2)
        # lnp += norm.logpdf(np.log(s), ) # TODO: put in hyper-parameters
        # TODO:
        pass

    return lnp

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))

def test_aic_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)
    aics = c.comparison.aic()
    assert len(aics) == 1
    assert aics[0] is None

def test_aic_fail_no_num_params():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_eff_data_points=1000)
    aics = c.comparison.aic()
    assert len(aics) == 1
    assert aics[0] is None

def test_aic_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)
    aics = c.comparison.aic()
    assert len(aics) == 1
    assert aics[0] == 0

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)

def test_aic_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)
    aics = c.comparison.aic()
    assert len(aics) == 2
    assert aics[0] == 0
    expected = 1 * 2 + (2.0 * 2 * 3 / (1000 - 2 - 1)) - (2.0 * 1 * 2 / (1000 - 1 - 1))
    assert np.isclose(aics[1], expected, atol=1e-3)

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

def test_bic_fail_no_num_params():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_eff_data_points=1000)
    bics = c.comparison.bic()
    assert len(bics) == 1
    assert bics[0] is None

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

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)

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)

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)

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

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)

def test_lpdf(x, loc, scale):
    aae(lpdf(x, loc, scale), norm.logpdf(x, loc, scale))

def test_lpdf_1d(x, loc, scale):
    aae(lpdf_1d(x, loc, scale), norm.logpdf(x, loc, scale))

def test_lpdf_3d(x, loc, scale):
    aae(lpdf_3d(x, loc, scale), norm.logpdf(x, loc, scale))

def test_lpdf_std(x):
    aae(lpdf_std(x), norm.logpdf(x))

def test_preds_ll(alpha, mu, gamma, err, num, w):
    current_impl = Lvm.preds_ll(alpha, mu, gamma, err, num, w)
    simple_impl = np.nansum(w * norm.logpdf(num, mu+gamma*alpha, err))
    simple_impl += np.sum(norm.logpdf(alpha))
    assert_approx_equal(current_impl, simple_impl)

def phoDurDist(x,dur_mean,proportion_std = 0.35):
    '''
    build gaussian distribution which has mean = dur_mean, std = dur_mean*proportion_std
    :param dur_mean: estimated from the centroid duration
    :param proportion_std:
    :return:
    '''
    prob = norm.logpdf(x,dur_mean,dur_mean*proportion_std)
    return prob

def _create_local_model_joint(self, targets, dataset):
        assert all(q in self.outputs for q in targets)
        assert dataset.shape[1] == len(targets)
        lookup = {
            'numerical': self._create_local_model_numerical,
            'categorical': self._create_local_model_categorical,
            'nominal': self._create_local_model_categorical,
        }
        models = {
            q: lookup[self.stattypes[self.outputs.index(q)]](q, dataset[:,i])
            for i, q in enumerate(targets)}
        simulate = lambda q, N=None: {c: models[c].simulate(N) for c in q}
        logpdf = lambda q: sum(models[c].logpdf(x) for c,x in q.iteritems())
        return LocalGpm(simulate, logpdf)

def _create_local_model_numerical(self, q, locality):
        assert q not in self.levels
        (mu, std) = (np.mean(locality), max(np.std(locality), .01))
        simulate = lambda N=None: self.rng.normal(mu, std, size=N)
        logpdf = lambda x: norm.logpdf(x, mu, std)
        return LocalGpm(simulate, logpdf)

def _create_local_model_categorical(self, q, locality):
        assert q in self.levels
        assert all(0 <= l < self.levels[q] for l in locality)
        counts = np.bincount(locality.astype(int), minlength=self.levels[q])
        p = counts / np.sum(counts, dtype=float)
        simulate = lambda N: self.rng.choice(self.levels[q], p=p, size=N)
        logpdf = lambda x: np.log(p[x])
        return LocalGpm(simulate, logpdf)

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)])

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]]))

def logpdf_marginal(self, z):
        return norm.logpdf(z, scale=1)

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))

def logpdf(self, rowid, targets, constraints=None, inputs=None):
        DistributionGpm.logpdf(self, rowid, targets, constraints, inputs)
        x = targets[self.outputs[0]]
        if not (self.l <= x <= self.h):
            return -float('inf')
        logpdf_unorm = NormalTrunc.calc_predictive_logp(
            x, self.mu, self.sigma, self.l, self.h)
        logcdf_norm = NormalTrunc.calc_log_normalizer(
            self.mu, self.sigma, self.l, self.h)
        return logpdf_unorm - logcdf_norm

def calc_predictive_logp(x, mu, sigma, l, h):
        return norm.logpdf(x, loc=mu, scale=sigma)

def __init__(self, K, Y, init=None, threshold=1e-9):

        N = np.shape(K)[0]
        f = np.zeros((N,1))
        converged = False
        k = 0
        innerC = 0

        for i in xrange(N):
            pdfDiff = norm.logpdf(f) - norm.logcdf(Y*f)
            W = np.exp(2*pdfDiff) + Y*f*np.exp(pdfDiff)
            Wsqrt = np.sqrt(W)
            Wdiag= np.diag(Wsqrt.flatten())

            B = np.identity(N) + np.dot(Wdiag, np.dot(K, Wdiag))
            grad = Y*np.exp(pdfDiff)
            b = W*f + grad
            interim = np.dot(Wdiag, np.dot(K, b))

            cgRes = Cg(B, interim, threshold=threshold)
            s1 = cgRes.result
            innerC = innerC + cgRes.iterations
            a = b - Wsqrt*s1

            if(converged):
                break
            f_prev = f
            f = np.dot(K, a)
            diff = f - f_prev
            if (np.dot(diff.T,diff).flatten() < threshold*N or innerC>15000):
                converged = True
            k = k+1

        self.result = f
        self.iterations = k + innerC

def __init__(self, K, Y, P, init=None, threshold=1e-9, precon=None):

        N = np.shape(K)[0]
        f = np.zeros((N,1))
        converged = False
        k = 0
        innerC = 0


        for i in xrange(N):
            pdfDiff = norm.logpdf(f) - norm.logcdf(Y*f)
            W = np.exp(2*pdfDiff) + Y*f*np.exp(pdfDiff)
            Wsqrt = np.sqrt(W)
            Wdiag= np.diag(Wsqrt.flatten())

            B = np.identity(N) + np.dot(Wdiag, np.dot(K, Wdiag))
            grad = Y*np.exp(pdfDiff)
            b = W*f + grad
            interim = np.dot(Wdiag, np.dot(K, b))

            pcgRes = RegularPcg(B, interim, None, threshold=threshold, preconInv=P.get_laplace_inversion(W,Wsqrt))
            s1 = pcgRes.result
            innerC = innerC + pcgRes.iterations
            a = b - Wsqrt*s1

            if(converged):
                break
            f_prev = f
            f = np.dot(K, a)
            diff = f - f_prev
            if (np.dot(diff.T,diff).flatten() < threshold*N or innerC>15000):
                converged = True
            k = k+1

        self.result = f
        self.iterations = k + innerC

def __init__(self, K, Y, init=None, threshold=1e-9):

        N = np.shape(K)[0]
        f = np.zeros((N,1))
        converged = False
        k = 0
        innerC = 0

        for i in xrange(N):
            pdfDiff = norm.logpdf(f) - norm.logcdf(Y*f)
            W = np.exp(2*pdfDiff) + Y*f*np.exp(pdfDiff)
            Wsqrt = np.sqrt(W)
            Wdiag= np.diag(Wsqrt.flatten())

            B = np.identity(N) + np.dot(Wdiag, np.dot(K, Wdiag))
            grad = Y*np.exp(pdfDiff)
            b = W*f + grad
            interim = np.dot(Wdiag, np.dot(K, b))

            cgRes = Cg(B, interim, threshold=threshold)
            s1 = cgRes.result
            innerC = innerC + cgRes.iterations
            a = b - Wsqrt*s1

            if(converged):
                break
            f_prev = f
            f = np.dot(K, a)
            diff = f - f_prev
            if (np.dot(diff.T,diff).flatten() < threshold*N or innerC>15000):
                converged = True
            k = k+1

        self.result = f
        self.iterations = k + innerC

def __init__(self, K, Y, P, init=None, threshold=1e-9, precon=None):

        N = np.shape(K)[0]
        f = np.zeros((N,1))
        converged = False
        k = 0
        innerC = 0


        for i in xrange(N):
            pdfDiff = norm.logpdf(f) - norm.logcdf(Y*f)
            W = np.exp(2*pdfDiff) + Y*f*np.exp(pdfDiff)
            Wsqrt = np.sqrt(W)
            Wdiag= np.diag(Wsqrt.flatten())

            B = np.identity(N) + np.dot(Wdiag, np.dot(K, Wdiag))
            grad = Y*np.exp(pdfDiff)
            b = W*f + grad
            interim = np.dot(Wdiag, np.dot(K, b))

            pcgRes = RegularPcg(B, interim, None, threshold=threshold, preconInv=P.get_laplace_inversion(W,Wsqrt))
            s1 = pcgRes.result
            innerC = innerC + pcgRes.iterations
            a = b - Wsqrt*s1

            if(converged):
                break
            f_prev = f
            f = np.dot(K, a)
            diff = f - f_prev
            if (np.dot(diff.T,diff).flatten() < threshold*N or innerC>15000):
                converged = True
            k = k+1

        self.result = f
        self.iterations = k + innerC