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