Python numpy.kron() Examples

The following are 30 code examples of numpy.kron(). 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 numpy , or try the search function .
Example #1
Source File: test_slinalg.py    From D-VAE with MIT License 6 votes vote down vote up
def test_perform(self):
        if not imported_scipy:
            raise SkipTest('kron tests need the scipy package to be installed')

        for shp0 in [(2,), (2, 3), (2, 3, 4), (2, 3, 4, 5)]:
            x = tensor.tensor(dtype='floatX',
                              broadcastable=(False,) * len(shp0))
            a = numpy.asarray(self.rng.rand(*shp0)).astype(config.floatX)
            for shp1 in [(6,), (6, 7), (6, 7, 8), (6, 7, 8, 9)]:
                if len(shp0) + len(shp1) == 2:
                    continue
                y = tensor.tensor(dtype='floatX',
                                  broadcastable=(False,) * len(shp1))
                f = function([x, y], kron(x, y))
                b = self.rng.rand(*shp1).astype(config.floatX)
                out = f(a, b)
                # Newer versions of scipy want 4 dimensions at least,
                # so we have to add a dimension to a and flatten the result.
                if len(shp0) + len(shp1) == 3:
                    scipy_val = scipy.linalg.kron(
                        a[numpy.newaxis, :], b).flatten()
                else:
                    scipy_val = scipy.linalg.kron(a, b)
                utt.assert_allclose(out, scipy_val) 
Example #2
Source File: test_gee.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def test_margins_gaussian(self):
        # Check marginal effects for a Gaussian GEE fit.  Marginal
        # effects and ordinary effects should be equal.

        n = 40
        np.random.seed(34234)
        exog = np.random.normal(size=(n, 3))
        exog[:, 0] = 1

        groups = np.kron(np.arange(n / 4), np.r_[1, 1, 1, 1])
        endog = exog[:, 1] + np.random.normal(size=n)

        model = sm.GEE(endog, exog, groups)
        result = model.fit(
            start_params=[-4.88085602e-04, 1.18501903, 4.78820100e-02])

        marg = result.get_margeff()

        assert_allclose(marg.margeff, result.params[1:])
        assert_allclose(marg.margeff_se, result.bse[1:])

        # smoke test
        marg.summary() 
Example #3
Source File: test_gee.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def test_constraint_covtype(self):
        # Test constraints with different cov types
        np.random.seed(6432)
        n = 200
        exog = np.random.normal(size=(n, 4))
        endog = exog[:, 0] + exog[:, 1] + exog[:, 2]
        endog += 3 * np.random.normal(size=n)
        group = np.kron(np.arange(n / 4), np.ones(4))
        L = np.array([[1., -1, 0, 0]])
        R = np.array([0., ])
        family = Gaussian()
        va = Independence()
        for cov_type in "robust", "naive", "bias_reduced":
            model = GEE(endog, exog, group, family=family,
                        cov_struct=va, constraint=(L, R))
            result = model.fit(cov_type=cov_type)
            result.standard_errors(cov_type=cov_type)
            assert_allclose(result.cov_robust.shape, np.r_[4, 4])
            assert_allclose(result.cov_naive.shape, np.r_[4, 4])
            if cov_type == "bias_reduced":
                assert_allclose(result.cov_robust_bc.shape, np.r_[4, 4]) 
Example #4
Source File: test_gee.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def test_linear_constrained(self):

        family = Gaussian()

        exog = np.random.normal(size=(300, 4))
        exog[:, 0] = 1
        endog = np.dot(exog, np.r_[1, 1, 0, 0.2]) +\
            np.random.normal(size=300)
        group = np.kron(np.arange(100), np.r_[1, 1, 1])

        vi = Independence()
        ve = Exchangeable()

        L = np.r_[[[0, 0, 0, 1]]]
        R = np.r_[0, ]

        for j, v in enumerate((vi, ve)):
            md = GEE(endog, exog, group, None, family, v,
                     constraint=(L, R))
            mdf = md.fit()
            assert_almost_equal(mdf.params[3], 0, decimal=10) 
Example #5
Source File: simulate_trotter_test.py    From OpenFermion-Cirq with Apache License 2.0 6 votes vote down vote up
def test_simulate_trotter_simulate_controlled(
        hamiltonian, time, initial_state, exact_state, order, n_steps,
        algorithm, result_fidelity):

    n_qubits = openfermion.count_qubits(hamiltonian)
    qubits = cirq.LineQubit.range(n_qubits)

    control = cirq.LineQubit(-1)
    zero = [1, 0]
    one = [0, 1]
    start_state = (numpy.kron(zero, initial_state)
                   + numpy.kron(one, initial_state)) / numpy.sqrt(2)

    circuit = cirq.Circuit(simulate_trotter(
        qubits, hamiltonian, time, n_steps, order, algorithm, control))

    final_state = circuit.final_wavefunction(start_state)
    correct_state = (numpy.kron(zero, initial_state)
                     + numpy.kron(one, exact_state)) / numpy.sqrt(2)
    assert fidelity(final_state, correct_state) > result_fidelity
    # Make sure the time wasn't too small
    assert fidelity(final_state, start_state) < 0.95 * result_fidelity 
Example #6
Source File: analysis.py    From OpenFermion-Cirq with Apache License 2.0 6 votes vote down vote up
def energy_from_opdm(opdm, constant, one_body_tensor, two_body_tensor):
    """
    Evaluate the energy of an opdm assuming the 2-RDM is opdm ^ opdm

    :param opdm: single spin-component of the full spin-orbital opdm.
    :param constant: constant shift to the Hamiltonian. Commonly this is the
                     nuclear repulsion energy.
    :param one_body_tensor: spatial one-body integrals
    :param two_body_tensor: spatial two-body integrals
    :return:
    """
    spin_opdm = np.kron(opdm, np.eye(2))
    spin_tpdm = 2 * wedge(spin_opdm, spin_opdm, (1, 1), (1, 1))
    molecular_hamiltonian = generate_hamiltonian(constant=constant,
                                                 one_body_integrals=one_body_tensor,
                                                 two_body_integrals=two_body_tensor)
    rdms = InteractionRDM(spin_opdm, spin_tpdm)
    return rdms.expectation(molecular_hamiltonian).real 
Example #7
Source File: mfopt.py    From OpenFermion-Cirq with Apache License 2.0 6 votes vote down vote up
def non_redundant_rotation_generators(
        rhf_objective: RestrictedHartreeFockObjective) -> List[FermionOperator]:  # testpragma: no cover
    # coverage: ignore
    """
    Generate the fermionic representation of all non-redundant rotation
    generators for restricted Hartree-fock

    :param rhf_objective: openfermioncirq.experiments.hfvqe.RestrictedHartreeFock object
    :return: list of fermionic generators.
    """
    rotation_generators = []
    for p in range(rhf_objective.nocc * rhf_objective.nvirt):
        grad_params = np.zeros(rhf_objective.nocc * rhf_objective.nvirt)
        grad_params[p] = 1
        kappa_spatial_orbital = rhf_params_to_matrix(grad_params,
                                                     len(rhf_objective.occ) + len(rhf_objective.virt),
                                                     rhf_objective.occ,
                                                     rhf_objective.virt)
        p0 = np.array([[1, 0], [0, 1]])
        kappa_spin_orbital = np.kron(kappa_spatial_orbital, p0)
        fermion_op = get_one_body_fermion_operator(kappa_spin_orbital)
        rotation_generators.append(fermion_op)
    return rotation_generators 
Example #8
Source File: test_gee_glm.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def setup_class(cls):

        vs = Independence()
        family = families.Gaussian()
        np.random.seed(987126)
        Y = np.random.normal(size=100)
        X1 = np.random.normal(size=100)
        X2 = np.random.normal(size=100)
        X3 = np.random.normal(size=100)
        groups = np.kron(np.arange(20), np.ones(5))

        D = pd.DataFrame({"Y": Y, "X1": X1, "X2": X2, "X3": X3})

        md = GEE.from_formula("Y ~ X1 + X2 + X3", groups, D,
                              family=family, cov_struct=vs)
        cls.result1 = md.fit()

        cls.result2 = GLM.from_formula("Y ~ X1 + X2 + X3", data=D).fit() 
Example #9
Source File: test_phreg.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def test_get_distribution(self):
        # Smoke test
        np.random.seed(34234)
        n = 200
        exog = np.random.normal(size=(n, 2))
        lin_pred = exog.sum(1)
        elin_pred = np.exp(-lin_pred)
        time = -elin_pred * np.log(np.random.uniform(size=n))
        status = np.ones(n)
        status[0:20] = 0
        strata = np.kron(range(5), np.ones(n // 5))

        mod = PHReg(time, exog, status=status, strata=strata)
        rslt = mod.fit()

        dist = rslt.get_distribution()

        fitted_means = dist.mean()
        true_means = elin_pred
        fitted_var = dist.var()
        fitted_sd = dist.std()
        sample = dist.rvs() 
Example #10
Source File: basistools.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def state_to_stdmx(state_vec):
    """
    Convert a state vector into a density matrix.

    Parameters
    ----------
    state_vec : list or tuple
       State vector in the standard (sigma-z) basis.

    Returns
    -------
    numpy.ndarray
        A density matrix of shape (d,d), corresponding to the pure state
        given by the length-`d` array, `state_vec`.
    """
    st_vec = state_vec.view(); st_vec.shape = (len(st_vec), 1)  # column vector
    dm_mx = _np.kron(_np.conjugate(_np.transpose(st_vec)), st_vec)
    return dm_mx  # density matrix in standard (sigma-z) basis 
Example #11
Source File: var_model.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def _var_acf(coefs, sig_u):
    """
    Compute autocovariance function ACF_y(h) for h=1,...,p

    Notes
    -----
    Lütkepohl (2005) p.29
    """
    p, k, k2 = coefs.shape
    assert(k == k2)

    A = util.comp_matrix(coefs)
    # construct VAR(1) noise covariance
    SigU = np.zeros((k*p, k*p))
    SigU[:k, :k] = sig_u

    # vec(ACF) = (I_(kp)^2 - kron(A, A))^-1 vec(Sigma_U)
    vecACF = scipy.linalg.solve(np.eye((k*p)**2) - np.kron(A, A), vec(SigU))

    acf = unvec(vecACF)
    acf = acf[:k].T.reshape((p, k, k))

    return acf 
Example #12
Source File: irf.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def _orth_cov(self):
        # Lutkepohl 3.7.8

        Ik = np.eye(self.neqs)
        PIk = np.kron(self.P.T, Ik)
        H = self.H

        covs = self._empty_covm(self.periods + 1)
        for i in range(self.periods + 1):
            if i == 0:
                apiece = 0
            else:
                Ci = np.dot(PIk, self.G[i-1])
                apiece = chain_dot(Ci, self.cov_a, Ci.T)

            Cibar = np.dot(np.kron(Ik, self.irfs[i]), H)
            bpiece = chain_dot(Cibar, self.cov_sig, Cibar.T) / self.T

            # Lutkepohl typo, cov_sig correct
            covs[i] = apiece + bpiece

        return covs 
Example #13
Source File: optools.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def unitary_to_process_mx(U):
    """
    Compute the super-operator which acts on (row)-vectorized
    density matrices from a unitary operator (matrix) U which
    acts on state vectors.  This super-operator is given by
    the tensor product of U and conjugate(U), i.e. kron(U,U.conj).

    Parameters
    ----------
    U : numpy array
        The unitary matrix which acts on state vectors.

    Returns
    -------
    numpy array
       The super-operator process matrix.
    """
    # U -> kron(U,Uc) since U rho U_dag -> kron(U,Uc)
    #  since AXB --row-vectorize--> kron(A,B.T)*vec(X)
    return _np.kron(U, _np.conjugate(U)) 
Example #14
Source File: vis_utils.py    From classification-of-encrypted-traffic with MIT License 6 votes vote down vote up
def plt_vector(x, colormap, num_headers):
    N = len(x)
    assert (N <= 16)
    len_x = 54
    len_y = num_headers
    # size = int(np.ceil(np.sqrt(len(x[0]))))
    length = len_y*len_x
    data = np.zeros((N, length), dtype=np.float64)
    data[:, :x.shape[1]] = x
    data = colormap(data / np.abs(data).max())
    # data = data.reshape([1, N, size, size, 3])
    data = data.reshape([1, N, len_y, len_x, 3])
    # data = np.pad(data, ((0, 0), (0, 0), (2, 2), (2, 2), (0, 0)), 'constant', constant_values=1)
    data = data.transpose([0, 2, 1, 3, 4]).reshape([1 * (len_y), N * (len_x), 3])
    return data
    # data = np.kron(data, np.ones([2, 2, 1])) # scales 
Example #15
Source File: vis_utils.py    From classification-of-encrypted-traffic with MIT License 6 votes vote down vote up
def graymap(x):

    x = x[..., np.newaxis]
    return np.concatenate([x, x, x], axis=-1)*0.5+0.5

# --------------------------------------
# Visualizing data
# --------------------------------------

# def visualize(x,colormap,name):
#
#     N = len(x)
#     assert(N <= 16)
#
#     x = colormap(x/np.abs(x).max())
#
#     # Create a mosaic and upsample
#     x = x.reshape([1, N, 29, 29, 3])
#     x = np.pad(x, ((0, 0), (0, 0), (2, 2), (2, 2), (0, 0)), 'constant', constant_values=1)
#     x = x.transpose([0, 2, 1, 3, 4]).reshape([1*33, N*33, 3])
#     x = np.kron(x, np.ones([2, 2, 1]))
#
#     PIL.Image.fromarray((x*255).astype('byte'), 'RGB').save(name) 
Example #16
Source File: test_kronecker.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up
def test_Kroneker(par):
    """Dot-test, inversion and comparison with np.kron for Kronecker operator
    """
    np.random.seed(10)
    G1 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    G2 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    x = np.ones(par['nx']**2) + par['imag']*np.ones(par['nx']**2)

    Kop = Kronecker(MatrixMult(G1, dtype=par['dtype']),
                    MatrixMult(G2, dtype=par['dtype']),
                    dtype=par['dtype'])
    assert dottest(Kop, par['ny']**2, par['nx']**2,
                   complexflag=0 if par['imag'] == 0 else 3)

    xlsqr = lsqr(Kop, Kop * x, damp=1e-20, iter_lim=300, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=2)

    # Comparison with numpy
    assert_array_almost_equal(np.kron(G1, G2), Kop * np.eye(par['nx']**2),
                              decimal=3) 
Example #17
Source File: b_model.py    From tenpy with GNU General Public License v3.0 6 votes vote down vote up
def init_H_bonds(self):
        """Initialize `H_bonds` hamiltonian.

        Called by __init__().
        """
        sx, sz, id = self.sigmax, self.sigmaz, self.id
        d = self.d
        nbonds = self.L - 1 if self.bc == 'finite' else self.L
        H_list = []
        for i in range(nbonds):
            gL = gR = 0.5 * self.g
            if self.bc == 'finite':
                if i == 0:
                    gL = self.g
                if i + 1 == self.L - 1:
                    gR = self.g
            H_bond = -self.J * np.kron(sx, sx) - gL * np.kron(sz, id) - gR * np.kron(id, sz)
            # H_bond has legs ``i, j, i*, j*``
            H_list.append(np.reshape(H_bond, [d, d, d, d]))
        self.H_bonds = H_list

    # (note: not required for TEBD) 
Example #18
Source File: tdvp_numpy.py    From tenpy with GNU General Public License v3.0 6 votes vote down vote up
def middle_bond_hamiltonian(Jx, Jz, hx, hz, L):
    """" Returns the spin operators sigma_x and sigma_z for L sites."""
    sx = np.array([[0., 1.], [1., 0.]])
    sz = np.array([[1., 0.], [0., -1.]])

    H_bond = Jx * np.kron(sx, sx) + Jz * np.kron(sz, sz)
    H_bond = H_bond + hx / 2 * np.kron(sx, np.eye(2)) + hx / 2 * np.kron(np.eye(2), sx)
    H_bond = H_bond + hz / 2 * np.kron(sz, np.eye(2)) + hz / 2 * np.kron(np.eye(2), sz)
    H_bond = H_bond.reshape(2, 2, 2, 2).transpose(0, 2, 1, 3).reshape(4, 4)  #i1 i2 i1' i2' -->
    U, s, V = np.linalg.svd(H_bond)

    M1 = np.dot(U, np.diag(s)).reshape(2, 2, 1, 4).transpose(2, 3, 0, 1)
    M2 = V.reshape(4, 1, 2, 2)
    M0 = np.tensordot(np.tensordot([1], [1], axes=0), np.eye(2), axes=0)
    W = []

    for i in range(L):
        if i == L / 2 - 1:
            W.append(M1)
        elif i == L / 2:
            W.append(M2)
        else:
            W.append(M0)
    return W 
Example #19
Source File: utils.py    From sparselandtools with MIT License 6 votes vote down vote up
def random_dictionary(n, K, normalized=True, seed=None):
    """
    Build a random dictionary matrix with K = n
    Args:
        n: square of signal dimension
        K: square of desired dictionary atoms
        normalized: If true, columns will be l2-normalized
        seed: Random seed

    Returns:
        Random dictionary
    """
    if seed:
        np.random.seed(seed)
    H = np.random.rand(n, K) * 255
    if normalized:
        for k in range(K):
            H[:, k] *= 1 / np.linalg.norm(H[:, k])
    return np.kron(H, H) 
Example #20
Source File: utils.py    From sparselandtools with MIT License 6 votes vote down vote up
def dictionary_from_transform(transform, n, K, normalized=True, inverse=True):
    """
    Builds a Dictionary matrix from a given transform

    Args:
        transform: A valid transform (e.g. Haar, DCT-II)
        n: number of rows transform dictionary
        K: number of columns transform dictionary
        normalized: If True, the columns will be l2-normalized
        inverse: Uses the inverse transform (as usually needed in applications)

    Returns:
        Dictionary build from the Kronecker-Delta of the transform applied to the identity.
    """
    H = np.zeros((K, n))
    for i in range(n):
        v = np.zeros(n)
        v[i] = 1.0
        H[:, i] = transform(v, sampling_factor=K)
    if inverse:
        H = H.T
    return np.kron(H.T, H.T) 
Example #21
Source File: vecm.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def cov_params_default(self):  # p.296 (7.2.21)
        # Sigma_co described on p. 287
        beta = self.beta
        if self.det_coef_coint.size > 0:
            beta = vstack((beta, self.det_coef_coint))
        dt = self.deterministic
        num_det = ("co" in dt) + ("lo" in dt)
        num_det += (self.seasons-1) if self.seasons else 0
        if self.exog is not None:
            num_det += self.exog.shape[1]
        b_id = scipy.linalg.block_diag(beta,
                                       np.identity(self.neqs * (self.k_ar-1) +
                                                   num_det))

        y_lag1 = self._y_lag1
        b_y = beta.T.dot(y_lag1)
        omega11 = b_y.dot(b_y.T)
        omega12 = b_y.dot(self._delta_x.T)
        omega21 = omega12.T
        omega22 = self._delta_x.dot(self._delta_x.T)
        omega = np.bmat([[omega11, omega12],
                         [omega21, omega22]]).A

        mat1 = b_id.dot(inv(omega)).dot(b_id.T)
        return np.kron(mat1, self.sigma_u) 
Example #22
Source File: irf.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def lr_effect_cov(self, orth=False):
        """
        Returns
        -------

        """
        lre = self.lr_effects
        Finfty = np.kron(np.tile(lre.T, self.lags), lre)
        Ik = np.eye(self.neqs)

        if orth:
            Binf = np.dot(np.kron(self.P.T, np.eye(self.neqs)), Finfty)
            Binfbar = np.dot(np.kron(Ik, lre), self.H)

            return (chain_dot(Binf, self.cov_a, Binf.T) +
                    chain_dot(Binfbar, self.cov_sig, Binfbar.T))
        else:
            return chain_dot(Finfty, self.cov_a, Finfty.T) 
Example #23
Source File: test_gee.py    From vnpy_crypto with MIT License 5 votes vote down vote up
def test_groups(self):
        # Test various group structures (nonconsecutive, different
        # group sizes, not ordered, string labels)

        n = 40
        x = np.random.normal(size=(n, 2))
        y = np.random.normal(size=n)

        # groups with unequal group sizes
        groups = np.kron(np.arange(n / 4), np.ones(4))
        groups[8:12] = 3
        groups[34:36] = 9

        model1 = GEE(y, x, groups=groups)
        result1 = model1.fit()

        # Unordered groups
        ix = np.random.permutation(n)
        y1 = y[ix]
        x1 = x[ix, :]
        groups1 = groups[ix]

        model2 = GEE(y1, x1, groups=groups1)
        result2 = model2.fit()

        assert_allclose(result1.params, result2.params)
        assert_allclose(result1.tvalues, result2.tvalues)

        # group labels are strings
        mp = {}
        import string
        for j, g in enumerate(set(groups)):
            mp[g] = string.ascii_letters[j:j + 4]
        groups2 = [mp[g] for g in groups]

        model3 = GEE(y, x, groups=groups2)
        result3 = model3.fit()

        assert_allclose(result1.params, result3.params)
        assert_allclose(result1.tvalues, result3.tvalues) 
Example #24
Source File: ar_model.py    From vnpy_crypto with MIT License 5 votes vote down vote up
def _presample_fit(self, params, start, p, end, y, predictedvalues):
        """
        Return the pre-sample predicted values using the Kalman Filter

        Notes
        -----
        See predict method for how to use start and p.
        """
        k = self.k_trend

        # build system matrices
        T_mat = KalmanFilter.T(params, p, k, p)
        R_mat = KalmanFilter.R(params, p, k, 0, p)

        # Initial State mean and variance
        alpha = np.zeros((p, 1))
        Q_0 = dot(inv(identity(p**2)-np.kron(T_mat, T_mat)),
                  dot(R_mat, R_mat.T).ravel('F'))

        Q_0 = Q_0.reshape(p, p, order='F')  # TODO: order might need to be p+k
        P = Q_0
        Z_mat = KalmanFilter.Z(p)
        for i in range(end):  # iterate p-1 times to fit presample
            v_mat = y[i] - dot(Z_mat, alpha)
            F_mat = dot(dot(Z_mat, P), Z_mat.T)
            Finv = 1./F_mat  # inv. always scalar
            K = dot(dot(dot(T_mat, P), Z_mat.T), Finv)
            # update state
            alpha = dot(T_mat, alpha) + dot(K, v_mat)
            L = T_mat - dot(K, Z_mat)
            P = dot(dot(T_mat, P), L.T) + dot(R_mat, R_mat.T)
            #P[0,0] += 1 # for MA part, R_mat.R_mat.T above
            if i >= start - 1:  # only record if we ask for it
                predictedvalues[i + 1 - start] = dot(Z_mat, alpha) 
Example #25
Source File: test_regression.py    From vnpy_crypto with MIT License 5 votes vote down vote up
def test_kron_matrix(self):
        # Ticket #71
        x = np.matrix('[1 0; 1 0]')
        assert_equal(type(np.kron(x, x)), type(x)) 
Example #26
Source File: test_phreg.py    From vnpy_crypto with MIT License 5 votes vote down vote up
def test_summary(self):
        # smoke test
        np.random.seed(34234)
        time = 50 * np.random.uniform(size=200)
        status = np.random.randint(0, 2, 200).astype(np.float64)
        exog = np.random.normal(size=(200,4))

        mod = PHReg(time, exog, status)
        rslt = mod.fit()
        smry = rslt.summary()

        strata = np.kron(np.arange(50), np.ones(4))
        mod = PHReg(time, exog, status, strata=strata)
        rslt = mod.fit()
        smry = rslt.summary()
        msg = "3 strata dropped for having no events"
        assert_(msg in str(smry))

        groups = np.kron(np.arange(25), np.ones(8))
        mod = PHReg(time, exog, status)
        rslt = mod.fit(groups=groups)
        smry = rslt.summary()

        entry = np.random.uniform(0.1, 0.8, 200) * time
        mod = PHReg(time, exog, status, entry=entry)
        rslt = mod.fit()
        smry = rslt.summary()
        msg = "200 observations have positive entry times"
        assert_(msg in str(smry)) 
Example #27
Source File: generalized_estimating_equations.py    From vnpy_crypto with MIT License 5 votes vote down vote up
def inverse(self, lpr):
        """
        Inverse of the multinomial logit transform, which gives the
        expected values of the data as a function of the linear
        predictors.

        Parameters
        ----------
        lpr : array-like (length must be divisible by `ncut`)
            The linear predictors

        Returns
        -------
        prob : array
            Probabilities, or expected values
        """

        expval = np.exp(lpr)

        denom = 1 + np.reshape(expval, (len(expval) // self.ncut,
                                        self.ncut)).sum(1)
        denom = np.kron(denom, np.ones(self.ncut, dtype=np.float64))

        prob = expval / denom

        return prob 
Example #28
Source File: slinalg.py    From D-VAE with MIT License 5 votes vote down vote up
def kron(a, b):
    """ Kronecker product.

    Same as scipy.linalg.kron(a, b).

    Parameters
    ----------
    a: array_like
    b: array_like

    Returns
    -------
    array_like with a.ndim + b.ndim - 2 dimensions

    Notes
    -----
    numpy.kron(a, b) != scipy.linalg.kron(a, b)!
    They don't have the same shape and order when
    a.ndim != b.ndim != 2.

    """
    a = tensor.as_tensor_variable(a)
    b = tensor.as_tensor_variable(b)
    if (a.ndim + b.ndim <= 2):
        raise TypeError('kron: inputs dimensions must sum to 3 or more. '
                        'You passed %d and %d.' % (a.ndim, b.ndim))
    o = tensor.outer(a, b)
    o = o.reshape(tensor.concatenate((a.shape, b.shape)),
                  a.ndim + b.ndim)
    shf = o.dimshuffle(0, 2, 1, * list(range(3, o.ndim)))
    if shf.ndim == 3:
        shf = o.dimshuffle(1, 0, 2)
        o = shf.flatten()
    else:
        o = shf.reshape((o.shape[0] * o.shape[2],
                         o.shape[1] * o.shape[3]) +
                        tuple(o.shape[i] for i in xrange(4, o.ndim)))
    return o 
Example #29
Source File: test_gee.py    From vnpy_crypto with MIT License 5 votes vote down vote up
def test_predict_exposure(self):

        n = 50
        X1 = np.random.normal(size=n)
        X2 = np.random.normal(size=n)
        groups = np.kron(np.arange(25), np.r_[1, 1])
        offset = np.random.uniform(1, 2, size=n)
        exposure = np.random.uniform(1, 2, size=n)
        Y = np.random.poisson(0.1 * (X1 + X2) + offset +
                              np.log(exposure), size=n)
        data = pd.DataFrame({"Y": Y, "X1": X1, "X2": X2, "groups": groups,
                             "offset": offset, "exposure": exposure})

        fml = "Y ~ X1 + X2"
        model = GEE.from_formula(fml, groups, data, family=Poisson(),
                                 offset="offset", exposure="exposure")
        result = model.fit()
        assert_equal(result.converged, True)

        pred1 = result.predict()
        pred2 = result.predict(offset=data["offset"])
        pred3 = result.predict(exposure=data["exposure"])
        pred4 = result.predict(
            offset=data["offset"], exposure=data["exposure"])
        pred5 = result.predict(exog=data[-10:],
                               offset=data["offset"][-10:],
                               exposure=data["exposure"][-10:])
        # without patsy
        pred6 = result.predict(exog=result.model.exog[-10:],
                               offset=data["offset"][-10:],
                               exposure=data["exposure"][-10:],
                               transform=False)
        assert_allclose(pred1, pred2)
        assert_allclose(pred1, pred3)
        assert_allclose(pred1, pred4)
        assert_allclose(pred1[-10:], pred5)
        assert_allclose(pred1[-10:], pred6) 
Example #30
Source File: test_slinalg.py    From D-VAE with MIT License 5 votes vote down vote up
def test_numpy_2d(self):
        for shp0 in [(2, 3)]:
            x = tensor.tensor(dtype='floatX',
                              broadcastable=(False,) * len(shp0))
            a = numpy.asarray(self.rng.rand(*shp0)).astype(config.floatX)
            for shp1 in [(6, 7)]:
                if len(shp0) + len(shp1) == 2:
                    continue
                y = tensor.tensor(dtype='floatX',
                                  broadcastable=(False,) * len(shp1))
                f = function([x, y], kron(x, y))
                b = self.rng.rand(*shp1).astype(config.floatX)
                out = f(a, b)
                assert numpy.allclose(out, numpy.kron(a, b))