update internal variables for sampling the distribution with the
        current covariance matrix C. This method is O(N^3), if C is not diagonal.
331913
def updateBD(self):
        
        # itereigenupdated is always up-to-date in the diagonal case
        # just double check here
        if self.itereigenupdated == self.countiter:
            return
        if self.opts['CMA_diagonal'] >= self.countiter:
            _print_warning("updateBD called in CMA_diagonal mode, " +
                           "this should be considered a bug", "updateBD",
                           iteration=self.countiter)
        # C has already positive updates, here come the additional negative updates
        if self.sp.neg.cmuexp:
            C_shrunken = (1 - self.sp.cmu - self.sp.c1)**(self.countiter - self.itereigenupdated)
            clip_fac = 0.60  # 0.9 is sufficient to prevent degeneration in small dimension
            if hasattr(self.opts['vv'], '__getitem__') and self.opts['vv'][0] == 'sweep_ccov_neg':
                clip_fac = 0.98
            if (self.countiter - self.itereigenupdated) * self.sp.neg.cmuexp * self.N \
                < clip_fac * C_shrunken:
                # pos.def. guarantied, because vectors are normalized
                self.C -= self.sp.neg.cmuexp * self._Yneg
            else:
                max_warns = 1
                try:
                    self._updateBD_warnings += 1
                except AttributeError:
                    self._updateBD_warnings = 1
                if self.opts['verbose'] > 1 and \
                                self._updateBD_warnings <= max_warns:
                    _print_warning('doing two additional eigen' +
                                   'decompositions to guarantee pos.def.',
                                   'updateBD', 'CMAEvolutionStrategy')
                    if self._updateBD_warnings == max_warns:
                        _print_warning('further warnings are surpressed',
                                       'updateBD')
                self.decompose_C()
                _tmp_inverse_root_C = dot(self.B / self.D, self.B.T)
                _tmp_inverse_root_C = (_tmp_inverse_root_C + _tmp_inverse_root_C.T) / 2
                Zneg = dot(dot(_tmp_inverse_root_C, self._Yneg), _tmp_inverse_root_C)
                eigvals, eigvecs = self.opts['CMA_eigenmethod'](Zneg)
                self.count_eigen += 1
                if max(eigvals) * self.sp.neg.cmuexp <= clip_fac:
                   self.C -= self.sp.neg.cmuexp * self._Yneg
                elif 1 < 3:
                    self.C -= (clip_fac / max(eigvals)) * self._Yneg
                    _print_warning(
                        'clipped learning rate for negative weights, ' +
                        'maximal eigenvalue = %f, maxeig * ccov = %f > %f'
                         % (max(eigvals), max(eigvals) * self.sp.neg.cmuexp, clip_fac),
                        iteration=self.countiter)
                    if 1 < 3: # let's check
                        eigvals, eigvecs = self.opts['CMA_eigenmethod'](self.C)
                        self.count_eigen += 1
                        print('new min eigenval = %e, old = %e'
                              % (min(eigvals), min(self.D)**2))
                        if min(eigvals) > 0:
                            print('new cond = %e, old = %e'
                                  % (max(eigvals) / min(eigvals),
                                     (max(self.D) / min(self.D))**2))
                else: # guaranties pos.def. unconditionally
                    _print_warning('exponential update for negative weights (internally more expensive)',
                                   iteration=self.countiter)
                    self.update_exponential(self._Yneg, -self.sp.neg.cmuexp)
                    # self.C = self.Ypos + Cs * Mh.expms(-self.sp.neg.cmuexp*Csi*self.Yneg*Csi) * Cs
            # Yneg = self.Yneg  # for temporary debugging, can be removed
            self._Yneg = np.zeros((self.N, self.N))

        if hasattr(self.opts['vv'], '__getitem__') and self.opts['vv'][0].startswith('sweep_ccov'):
            self.opts['CMA_const_trace'] = True
        if self.opts['CMA_const_trace'] in (True, 1, 2):  # normalize trace of C
            if self.opts['CMA_const_trace'] == 2:
                s = np.exp(2 * np.mean(np.log(self.D)))  # or geom average of dC?
            else:
                s = np.mean(np.diag(self.C))
            self.C /= s

        dC = np.diag(self.C)
        if max(dC) / min(dC) > 1e8:
            # allows for much larger condition numbers, if axis-parallel
            self.sigma_vec *= np.diag(self.C)**0.5
            self.C = self.correlation_matrix()
            _print_warning('condition in coordinate system exceeded 1e8' +
                ', rescaled to 1')

        # self.C = np.triu(self.C) + np.triu(self.C,1).T  # should work as well
        # self.D, self.B = eigh(self.C) # hermitian, ie symmetric C is assumed

        self.decompose_C()
        # assert(sum(self.D-DD) < 1e-6)
        # assert(sum(sum(np.dot(BB, BB.T)-np.eye(self.N))) < 1e-6)
        # assert(sum(sum(np.dot(BB * DD, BB.T) - self.C)) < 1e-6)
        # assert(all(self.B[self.countiter % self.N] == self.B[self.countiter % self.N,:]))

        # qqqqqqqqqq
        # is O(N^3)
        # assert(sum(abs(self.C - np.dot(self.D * self.B,  self.B.T))) < N**2*1e-11)

        if 1 < 3 and max(self.D) / min(self.D) > 1e6 and self.gp.isidentity:
            # TODO: allow to do this again
            # dmean_prev = dot(self.B, (1. / self.D) * dot(self.B.T, (self.mean - 0*self.mean_old) / self.sigma_vec))
            self.gp._tf_matrix = (self.sigma_vec * dot(self.B * self.D, self.B.T).T).T
            self.gp._tf_matrix_inv = (dot(self.B / self.D, self.B.T).T / self.sigma_vec).T
            self.gp.tf_pheno = lambda x: dot(self.gp._tf_matrix, x)
            self.gp.tf_geno = lambda x: dot(self.gp._tf_matrix_inv, x)  # not really necessary
            self.gp.isidentity = False
            assert self.mean is not self.mean_old
            self.mean = self.gp.geno(self.mean)  # same as tf_geno
            self.mean_old = self.gp.geno(self.mean_old)  # not needed?
            self.pc = self.gp.geno(self.pc)
            self.D[:] = 1.0
            self.B = np.eye(self.N)
            self.C = np.eye(self.N)
            self.dC[:] = 1.0
            self.sigma_vec = 1
            # dmean_now = dot(self.B, (1. / self.D) * dot(self.B.T, (self.mean - 0*self.mean_old) / self.sigma_vec))
            # assert Mh.vequals_approximately(dmean_now, dmean_prev)
            _print_warning('\n geno-pheno transformation introduced based on current C,\n injected solutions become "invalid" in this iteration',
                           'updateBD', 'CMAEvolutionStrategy', self.countiter)

        self.itereigenupdated = self.countiter
495777
def _update_estimate_and_sampler(self, ell, ell_hat, weight, extra_info,
                                     **kwargs):
        """"""
        self._TP += ell_hat * ell * weight
        self._FP += ell_hat * (1 - ell) * weight
        self._FN += (1 - ell_hat) * ell * weight
        self._TN += (1 - ell_hat) * (1 - ell) * weight

        self._estimate[self.t_] = \
                self._F_measure(self.alpha, self._TP, self._FP, self._FN)
697893
def sample_W(self, nSamples, raw_samples=False):
        
        assert isinstance(self.kern, kern.Linear)
        from ..util.linalg import pdinv
        N, D = self.Y.shape
        Q = self.X.shape[1]
        noise_var = self.likelihood.variance.values
        
        # Draw samples for X
        Xs = np.random.randn(*((nSamples,)+self.X.shape))*np.sqrt(self.X.variance.values)+self.X.mean.values
        b = np.random.rand(*((nSamples,)+self.X.shape))
        Xs[b>self.X.gamma.values] = 0
        
        invcov = (Xs[:,:,:,None]*Xs[:,:,None,:]).sum(1)/noise_var+np.eye(Q)
        cov = np.array([pdinv(invcov[s_idx])[0] for s_idx in range(invcov.shape[0])])
        Ws = np.empty((nSamples, Q, D))
        tmp = (np.transpose(Xs, (0,2,1)).reshape(nSamples*Q,N).dot(self.Y)).reshape(nSamples,Q,D)
        mean = (cov[:,:,:,None]*tmp[:,None,:,:]).sum(2)/noise_var
        zeros = np.zeros((Q,))
        for s_idx in range(Xs.shape[0]):
            Ws[s_idx] = (np.random.multivariate_normal(mean=zeros,cov=cov[s_idx],size=(D,))).T+mean[s_idx]
        
        if raw_samples:
            return Ws
        else:
            return Ws.mean(0), Ws.std(0)
1281247
function cov(x) {
	// Useful variables
	let N = x.length;

	// Step 1: Find expected value for each variable (columns are variables, rows are samples)
	let E = [];
	for (let c = 0; c < x[0].length; c++) {
		let u = 0;
		for (let r = 0; r < N; r++) {
			u += x[r][c];
		}
		E.push(u / N);
	}

	// Step 2: Center each variable at 0
	let centered = [];
	for (let r = 0; r < N; r++) {
		centered.push([]);
		for (let c = 0; c < x[0].length; c++) {
			centered[r].push(x[r][c] - E[c]);
		}
	}

	// Step 3: Calculate covariance between each possible combination of variables
	let S = [];

	for (let i = 0; i < x[0].length; i++) {
		S.push([]);
	}

	for (let c1 = 0; c1 < x[0].length; c1++) {
		// Doing c2 = c1; because, for example, cov(1,2) = cov(2, 1)
		// so no need to recalculate
		for (let c2 = c1; c2 < x[0].length; c2++) {
			// Solve cov(c1, c2) = average of elementwise products of each variable sample
			let cov = 0;
			for (let r = 0; r < N; r++) {
				cov += centered[r][c1] * centered[r][c2];
			}
			cov /= (N - 1); // N-1 for sample covariance, N for population. In this case, using sample covariance
			S[c1][c2] = cov;
			S[c2][c1] = cov; // Matrix is symmetric
		}
	}

	return S;
}
12453
def _update_cov_model(self, strata_to_update='all'):
        
        if strata_to_update == 'all':
            strata_to_update = self.strata.indices_
        #: Otherwise assume strata_to_update is valid (no duplicates etc.)

        #: Update covariance matrices
        #: We usually update only one stratum at a time, so for loop is ok
        n_sampled = np.clip(self.strata._n_sampled, 2, np.inf) #: adding 2 avoids undef. cov
        factor = n_sampled/(n_sampled - 1)
        for k in strata_to_update:
            TP = self._BB_TP.theta_[k]
            PP = self._BB_PP.theta_[k]
            P = self._BB_P.theta_[k]

            self.cov_model_[k,0,0] = factor[k] * TP * (1 - TP)
            self.cov_model_[k,0,1] = factor[k] * TP * (1 - PP)
            self.cov_model_[k,0,2] = factor[k] * TP * (1 - P)
            self.cov_model_[k,1,1] = factor[k] * PP * (1 - PP)
            self.cov_model_[k,1,2] = factor[k] * (TP - PP * P)
            self.cov_model_[k,2,2] = factor[k] * P * (1 - P)
            self.cov_model_[k,1,0] = self.cov_model_[k,0,1]
            self.cov_model_[k,2,0] = self.cov_model_[k,0,2]
            self.cov_model_[k,2,1] = self.cov_model_[k,1,2]
838795
def propose(self):
        

        # Locally store size of matrix
        dims = self.stochastic.value.shape

        # Add normal deviate to value and symmetrize
        dev = rnormal(
            0,
            self.adaptive_scale_factor *
            self.proposal_sd,
            size=dims)
        symmetrize(dev)

        # Replace
        self.stochastic.value = dev + self.stochastic.value
515222
def _forward_step(self, Y, T_X, Var_X, Var_dX, rho_e, sigma2_e, weight):
        
        # We currently only implement diagonal form
        # of covariance matrix for Var_X, Var_dX and T_X, which means
        # each dimension of X is independent and their innovation noise
        # are also independent. Note that log_p_data takes this assumption.
        if Var_X.ndim == 1:
            inv_Var_X = np.diag(1 / Var_X)
            half_log_det_Var_X = np.sum(np.log(Var_X)) / 2.0
            Var_X = np.diag(Var_X)
            # the marginal variance of X
        else:
            half_log_det_Var_X = self._half_log_det(Var_X)
            inv_Var_X = np.linalg.inv(Var_X)
        if Var_dX.ndim == 1:
            inv_Var_dX = np.diag(1 / Var_dX)
            half_log_det_Var_dX = np.sum(np.log(Var_dX)) / 2.0
            Var_dX = np.diag(Var_dX)
            # the marginal variance of Delta X (the change of X from
            # previous time point)
        else:
            inv_Var_dX = np.linalg.inv(Var_dX)
            half_log_det_Var_dX = self._half_log_det(Var_dX)
        if T_X.ndim == 1:
            T_X = np.diag(T_X)
            # Transfer function of X: the expected mean of X at t+1
            # time point is T_x * X
        [n_T, n_V] = np.shape(Y)
        # numbers of time points and voxels
        mu = [None] * n_T
        # posterior mean of X, conditioned on all data up till the current
        # time point
        Gamma_inv = [None] * n_T
        # inverse of poterior Gamma.
        mu_Gamma_inv = [None] * n_T
        # mu * inv(Gamma)
        log_p_data = - np.log(np.pi * 2) * (n_T * n_V) / 2 \
            - half_log_det_Var_X - np.sum(np.log(sigma2_e)) * n_T / 2.0\
            + np.sum(np.log(1 - rho_e**2)) / 2.0 - half_log_det_Var_dX \
            * (n_T - 1)
        # This is the term to be incremented by c_n at each time step.
        # We first add all the fixed terms to it.

        # The following are a few fixed terms.
        Lambda_0 = np.dot(T_X, np.dot(inv_Var_dX, T_X.T)) \
            + np.dot(weight * rho_e**2 / sigma2_e, weight.T)
        H = np.dot(inv_Var_dX, T_X.T) + np.dot(weight * rho_e / sigma2_e,
                                               weight.T)
        Lambda_1 = inv_Var_dX + np.dot(weight / sigma2_e, weight.T)

        Gamma_inv[0] = inv_Var_X + np.dot(
            weight * (1 - rho_e**2) / sigma2_e, weight.T)
        # We might not need this and only use linalg.solve for related terms.
        mu_Gamma_inv[0] = np.dot(
            Y[0, :] * (1 - rho_e**2) / sigma2_e, weight.T)
        mu[0] = np.linalg.solve(Gamma_inv[0], mu_Gamma_inv[0])
        log_p_data -= 0.5 * np.sum(Y[0, :]**2 * (1 - rho_e**2) / sigma2_e)
        # This is the term added for the first time point.

        deltaY = Y[1:, :] - rho_e * Y[:-1, :]
        deltaY_sigma2inv_rho_weightT = np.dot(
            deltaY / sigma2_e * rho_e, weight.T)
        for t in np.arange(1, n_T):
            Gamma_tilde_inv = Lambda_0 + Gamma_inv[t - 1]
            tmp = np.linalg.solve(Gamma_tilde_inv, H.T)
            Gamma_inv[t] = Lambda_1 - np.dot(H, tmp)
            mu_Gamma_inv[t] = np.dot(deltaY[t - 1, :] / sigma2_e, weight.T) \
                + np.dot(mu_Gamma_inv[t - 1]
                         - deltaY_sigma2inv_rho_weightT[t - 1, :], tmp)
            mu[t] = np.linalg.solve(Gamma_inv[t], mu_Gamma_inv[t])
            tmp2 = mu_Gamma_inv[t - 1] - deltaY_sigma2inv_rho_weightT[t - 1, :]
            log_p_data += -self._half_log_det(Gamma_tilde_inv) \
                + np.dot(tmp2, np.linalg.solve(Gamma_tilde_inv, tmp2)) / 2.0
        log_p_data += -self._half_log_det(Gamma_inv[-1]) \
            + np.dot(mu_Gamma_inv[-1], mu[-1]) / 2.0 \
            - np.sum(deltaY**2 / sigma2_e) / 2.0
        return mu, mu_Gamma_inv, Gamma_inv, log_p_data, Lambda_0, \
            Lambda_1, H, deltaY, deltaY_sigma2inv_rho_weightT
513363
def update_gradients_full(self, dL_dK, X, X2=None):
        """"""
        dvar, dw, db = self._comp_grads(dL_dK, X, X2)[:3]
        self.variance.gradient = dvar
        self.weight_variance.gradient = dw
        self.bias_variance.gradient = db
327156
def _refit(self):
        
        # temporary variables for manipulation
        mu_diff = np.empty_like(self.f.components[0].mu)
        sigma   = np.empty_like(self.f.components[0].sigma)
        mean    = np.empty_like(mu_diff)
        cov     = np.empty_like(sigma)

        for j, c in enumerate(self.g.components):
            # stop if inv_map is empty for j-th comp.
            if not self.inv_map[j]:
                self.g.weights[j] = 0.
                continue

            # (re-)initialize new mean/cov to zero
            mean[:] = 0.0
            cov[:] = 0.0

            # compute total weight and mean
            self.g.weights[j] = self.f.weights[self.inv_map[j]].sum()
            for i in self.inv_map[j]:
                mean += self.f.weights[i] * self.f.components[i].mu

            # rescale by total weight
            mean /= self.g.weights[j]

            # update covariance
            for i in self.inv_map[j]:
                # mu_diff = mu'_j - mu_i
                mu_diff[:] = mean
                mu_diff -= self.f.components[i].mu

                # sigma = (mu'_j - mu_i) (mu'_j - mu_i)^T
                sigma[:] = np.outer(mu_diff, mu_diff)

                # sigma += sigma_i
                sigma += self.f.components[i].sigma

                # multiply with alpha_i
                sigma *= self.f.weights[i]

                # sigma_j += alpha_i * (sigma_i + (mu'_j - mu_i) (mu'_j - mu_i)^T
                cov += sigma

            # 1 / beta_j
            cov /= self.g.weights[j]

            # update the Mixture
            c.update(mean, cov)
1359199
public static Matrix calculateCovarianceMatrixOfCenteredData(Matrix data){
		Assume.assume(data.getRowDimension()>1, "Number of data samples is "+
				data.getRowDimension()+", but it has to be >1 to compute "+
				"covariances");
		int dimsNo = data.getColumnDimension();
		int samplesNo = data.getRowDimension();
		Matrix m = new Matrix(dimsNo, dimsNo);
		for(int r = 0; r < dimsNo; r++)
			for(int c = r; c < dimsNo; c++){
				double sum = 0;
				for(int i = 0; i < samplesNo; i++)
					 sum += data.get(i, r)*data.get(i, c);
				m.set(r, c, sum/(samplesNo-1));
			}
		for(int r = 0; r < dimsNo; r++)
			for(int c = 0; c < r; c++) m.set(r, c, m.get(c, r));
		return m;		
	}
}
939687
def G(self):
        
        if self.c is None:
            raise RuntimeError('Inverse cross spectral density requires '
                               'invertible noise covariance matrix c.')
        A = self.A()
        # TODO: can we do that more efficiently?
        G = np.einsum('ji..., jk... ->ik...', A.conj(), self.Cinv())
        G = np.einsum('ij..., jk... ->ik...', G, A)
        return G
1149114
def _update_inernal(self, varp_list):
        """"""
        # The probability for the binary variable for the same latent dimension of any of the models is on.
        self._b_prob_all = 1.-param_to_array(varp_list[0].gamma_group)
        [np.multiply(self._b_prob_all, 1.-vp.gamma_group, self._b_prob_all) for vp in varp_list[1:]]
526360
def updateModel(self, X_all, Y_all, X_new, Y_new):
        
        if self.model is None:
            self._create_model(X_all, Y_all)
        else:
            self.model.set_XY(X_all, Y_all)

        # WARNING: Even if self.max_iters=0, the hyperparameters are bit modified...
        if self.max_iters > 0:
            # --- update the model maximizing the marginal likelihood.
            if self.optimize_restarts==1:
                self.model.optimize(optimizer=self.optimizer, max_iters = self.max_iters, messages=False, ipython_notebook=False)
            else:
                self.model.optimize_restarts(num_restarts=self.optimize_restarts, optimizer=self.optimizer, max_iters = self.max_iters, verbose=self.verbose)
269460
def shrink(self):
        
        x = np.nan_to_num(self.X.values)

        # de-mean returns
        t, n = np.shape(x)
        meanx = x.mean(axis=0)
        x = x - np.tile(meanx, (t, 1))
        xmkt = x.mean(axis=1).reshape(t, 1)

        # compute sample covariance matrix
        sample = np.cov(np.append(x, xmkt, axis=1), rowvar=False) * (t - 1) / t
        covmkt = sample[0:n, n].reshape(n, 1)
        varmkt = sample[n, n]
        sample = sample[:n, :n]
        prior = np.dot(covmkt, covmkt.T) / varmkt
        prior[np.eye(n) == 1] = np.diag(sample)

        # compute shrinkage parameters
        if self.delta is None:
            c = np.linalg.norm(sample - prior, "fro") ** 2
            y = x ** 2
            p = 1 / t * np.sum(np.dot(y.T, y)) - np.sum(sample ** 2)
            # r is divided into diagonal
            # and off-diagonal terms, and the off-diagonal term
            # is itself divided into smaller terms
            rdiag = 1 / t * np.sum(y ** 2) - sum(np.diag(sample) ** 2)
            z = x * np.tile(xmkt, (n,))
            v1 = 1 / t * np.dot(y.T, z) - np.tile(covmkt, (n,)) * sample
            roff1 = (
                np.sum(v1 * np.tile(covmkt, (n,)).T) / varmkt
                - np.sum(np.diag(v1) * covmkt.T) / varmkt
            )
            v3 = 1 / t * np.dot(z.T, z) - varmkt * sample
            roff3 = (
                np.sum(v3 * np.dot(covmkt, covmkt.T)) / varmkt ** 2
                - np.sum(np.diag(v3).reshape(-1, 1) * covmkt ** 2) / varmkt ** 2
            )
            roff = 2 * roff1 - roff3
            r = rdiag + roff

            # compute shrinkage constant
            k = (p - r) / c
            shrinkage = max(0, min(1, k / t))
            self.delta = shrinkage
        else:
            # use specified constant
            shrinkage = self.delta

        # compute the estimator
        sigma = shrinkage * prior + (1 - shrinkage) * sample
        return self.format_and_annualise(sigma)
1295057
def resample(self):
        

        if self.just_resampled:
            warnings.warn(
                "Resampling without additional data; this may not perform as "
                "desired.",
                ResamplerWarning
            )

        # Record that we have performed a resampling step.
        self._just_resampled = True
        self._resample_count += 1

        # If we're tracking divergences, make a copy of the weights and
        # locations.
        if self._resampling_divergences is not None:
            old_locs = self.particle_locations.copy()
            old_weights = self.particle_weights.copy()

        # Record the previous mean, cov if needed.
        if self._debug_resampling:
            old_mean = self.est_mean()
            old_cov = self.est_covariance_mtx()

        # Find the new particle locations according to the chosen resampling
        # algorithm.
        # We pass the model so that the resampler can check for validity of
        # newly placed particles.
        # FIXME This feels fishy. If we update particles elsewwhere
        new_distribution = self.resampler(self.model, self)
        self.particle_weights = new_distribution.particle_weights
        self.particle_locations = new_distribution.particle_locations

        # Possibly canonicalize, if we've been asked to do so.
        if self._canonicalize:
            self.particle_locations[:, :] = self.model.canonicalize(self.particle_locations)

        # Instruct the model to clear its cache, demoting any errors to
        # warnings.
        try:
            self.model.clear_cache()
        except Exception as e:
            warnings.warn("Exception raised when clearing model cache: {}. Ignoring.".format(e))

        # Possibly track the new divergence.
        if self._resampling_divergences is not None:
            self._resampling_divergences.append(
                self._kl_divergence(old_locs, old_weights)
            )

        # Report current and previous mean, cov.
        if self._debug_resampling:
            new_mean = self.est_mean()
            new_cov = self.est_covariance_mtx()
            logger.debug("Resampling changed mean by {}. Norm change in cov: {}.".format(
                old_mean - new_mean,
                np.linalg.norm(new_cov - old_cov)
            ))
172824
def update(self):
        '''
        
        '''
        observed_arr = self.noise_sampler.generate()
        inferenced_arr = self.inference(observed_arr)

        error_arr = self.__convolutional_auto_encoder.computable_loss.compute_loss(
            observed_arr,
            inferenced_arr
        )

        delta_arr = self.__convolutional_auto_encoder.computable_loss.compute_delta(
            observed_arr,
            inferenced_arr
        )

        delta_arr = self.__convolutional_auto_encoder.back_propagation(delta_arr)
        self.__convolutional_auto_encoder.optimize(self.__learning_rate, 1)

        return error_arr
1244175
def sample_niw(mu,lmbda,kappa,nu):
    '''
    
    '''
    # code is based on Matlab's method
    # reference: p. 87 in Gelman's Bayesian Data Analysis
    assert nu > lmbda.shape[0] and kappa > 0

    # first sample Sigma ~ IW(lmbda,nu)
    lmbda = sample_invwishart(lmbda,nu)
    # then sample mu | Lambda ~ N(mu, Lambda/kappa)
    mu = np.random.multivariate_normal(mu,lmbda / kappa)

    return mu, lmbda