Get the log2 bayes factor of the fit for each modality
338211
def single_feature_fit(self, feature):
        """"""
        if np.isfinite(feature).sum() == 0:
            series = pd.Series(index=MODALITY_ORDER)
        else:
            logbf_one_param = pd.Series(
                {k: v.logsumexp_logliks(feature) for
                 k, v in self.one_param_models.items()})

            # Check if none of the previous features fit
            if (logbf_one_param <= self.logbf_thresh).all():
                logbf_two_param = pd.Series(
                    {k: v.logsumexp_logliks(feature)
                     for k, v in self.two_param_models.items()})
                series = pd.concat([logbf_one_param, logbf_two_param])
                series[NULL_MODEL] = self.logbf_thresh
            else:
                series = logbf_one_param
        series.index.name = 'Modality'
        series.name = self.score_name
        return series
1184175
def do_rmfit(self, arg):
        """"""
        if arg in self.curargs["fits"]:
            del self.curargs["fits"][arg]        
        #We also need to remove the variable entry if it exists.
        if "timing" in arg:
            fitvar = "{}|fit".format(arg)
        else:
            fitvar = "{}.fit".format(arg)
        if fitvar in self.curargs["dependents"]:
            self.curargs["dependents"].remove(fitvar)
1180004
def get_logreg_coefs(self, category, clf=LogisticRegression()):
        ''' 
        '''
        self._fit_tfidf_model(category, clf)
        return clf.coef_[0]
1307427
def fmt_chi2(f):
    
    if hasattr(f, 'logGBF'):
        fmt = "chi2/dof = %.2g [%d]    Q = %.2g    log(GBF) = %.5g"
        chi2_dof = f.chi2 / f.dof if f.dof != 0 else 0
        return fmt % (chi2_dof, f.dof, f.Q, f.logGBF)
    else:
        fmt = "chi2/dof = %.2g [%d]    Q = %.2g"
        chi2_dof = f.chi2 / f.dof if f.dof != 0 else 0
        return fmt % (chi2_dof, f.dof, f.Q)
247732
def getCandScoresMapBruteForce(self, profile):
        

        wmg = profile.getWmg(True)
        m = len(wmg.keys())
        cands = range(m)
        V = self.createBinaryRelation(m)
        gains = dict()
        for cand in wmg.keys():
            gains[cand] = 0
        graphs = itertools.product(range(2), repeat=m*(m-1)/2)
        for comb in graphs:
            prob = 1
            i = 0
            for a, b in itertools.combinations(cands,2):
                V[a][b] = comb[i]
                V[b][a] = 1-comb[i]
                if comb[i] > 0:
                    prob *= 1/(1+self.phi ** float(wmg[a+1][b+1]))
                else:
                    prob *= 1/(1+self.phi ** float(wmg[b+1][a+1]))
                i += 1
                if i >= m*(m-1)/2:
                    break
            for cand in wmg.keys():
                gains[cand] += self.utilityFunction.getUtility([cand], V)*prob
        return gains
168856
def fit(self, X, y=None):
        
        if self.normalize:
            X = normalize(X)

        self._check_force_weights()
        random_state = check_random_state(self.random_state)
        X = self._check_fit_data(X)

        (
            self.cluster_centers_,
            self.labels_,
            self.inertia_,
            self.weights_,
            self.concentrations_,
            self.posterior_,
        ) = movMF(
            X,
            self.n_clusters,
            posterior_type=self.posterior_type,
            force_weights=self.force_weights,
            n_init=self.n_init,
            n_jobs=self.n_jobs,
            max_iter=self.max_iter,
            verbose=self.verbose,
            init=self.init,
            random_state=random_state,
            tol=self.tol,
            copy_x=self.copy_x,
        )

        return self
710441
def get_kl_loss(self, means, log_vars, means_p=None, log_vars_p=None):
    """"""
    kl_loss = 0.0
    if means_p is None:
      means_p = tf.unstack(tf.zeros_like(means))
    if log_vars_p is None:
      log_vars_p = tf.unstack(tf.zeros_like(log_vars))
    enumerated_inputs = enumerate(zip(means, log_vars, means_p, log_vars_p))
    if self.is_training and self.hparams.stochastic_model:
      for i, (mean, log_var, mean_p, log_var_p) in enumerated_inputs:
        kl_loss += common_layers.kl_divergence(mean, log_var, mean_p, log_var_p)
        tf.summary.histogram("posterior_mean_%d" % i, mean)
        tf.summary.histogram("posterior_log_var_%d" % i, log_var)
        tf.summary.histogram("prior_mean_%d" % i, mean_p)
        tf.summary.histogram("prior_log_var_%d" % i, log_var_p)
      tf.summary.scalar("kl_raw", tf.reduce_mean(kl_loss))

    beta = self.get_beta(kl_loss)
    # information capacity from "Understanding disentangling in beta-VAE"
    if self.hparams.information_capacity > 0.0:
      kl_loss = tf.abs(kl_loss - self.hparams.information_capacity)
    return beta * kl_loss
331430
def _getScalesPairwise(self,verbose=False, initDiagonal=False):
        
        var = sp.zeros((self.P,2))

        if initDiagonal:
            #1. fit single trait model
            if verbose:
                print('.. fit single-trait model for initialization')
            vc = VarianceDecomposition(self.Y[:,0:1])
            for term_i in range(self.n_randEffs):
                if term_i==self.noisPos:
                    vc.addRandomEffect(is_noise=True)
                else:
                    K = self.vd.getTerm(term_i).getK()
                    vc.addRandomEffect(K=K)
            scales0 = sp.sqrt(0.5)*sp.ones(2)

            for p in range(self.P):
                if verbose: print(('   .. trait %d' % p))
                vc.setY(self.Y[:,p:p+1])
                conv = vc.optimize(scales0=scales0)
                if not conv:
                    print('warning initialization not converged')
                var[p,:] = vc.getVarianceComps()[0,:]

        elif fastlmm_present:
            if verbose:
                print('.. fit single-trait model for initialization (using fastlmm)')
            for p in range(self.P):
                if verbose: print(('   .. trait %d' % p))
                covariates = None
                for term_i in range(self.n_randEffs):
                    if term_i==self.noisPos:
                        pass
                    else:
                        K = self.vd.getTerm(term_i).getK()
                varY = sp.var(self.Y[:,p:p+1])
                lmm = fastLMM(X=covariates, Y=self.Y[:,p:p+1], G=None, K=K)
                opt = lmm.findH2(nGridH2=100)
                h2 = opt['h2']
                var[p,:] = h2 * varY
                var[p,self.noisPos] = (1.0-h2) * varY
                #import ipdb;ipdb.set_trace()
        else:
            if verbose:
                print('.. random initialization of diagonal')
            var = sp.random.randn(var.shape[0],var.shape[1])
            var = var*var + 0.001
        #2. fit pairwise model
        if verbose:
            print('.. fit pairwise model for initialization')
        vc = VarianceDecomposition(self.Y[:,0:2])
        for term_i in range(self.n_randEffs):
            if term_i==self.noisPos:
                vc.addRandomEffect(is_noise=True,trait_covar_type='freeform')
            else:
                K = self.vd.getTerm(term_i).getK()
                vc.addRandomEffect(K=K,trait_covar_type='freeform')
        rho_g = sp.ones((self.P,self.P))
        rho_n = sp.ones((self.P,self.P))
        for p1 in range(self.P):
            for p2 in range(p1):
                if verbose:
                    print(('   .. fit pair (%d,%d)'%(p1,p2)))
                vc.setY(self.Y[:,[p1,p2]])
                scales0 = sp.sqrt(sp.array([var[p1,0],1e-4,var[p2,0],1e-4,var[p1,1],1e-4,var[p2,1],1e-4]))
                conv = vc.optimize(scales0=scales0)
                if not conv:
                    print('warning initialization not converged')
                Cg = vc.getTraitCovar(0)
                Cn = vc.getTraitCovar(1)
                rho_g[p1,p2] = Cg[0,1]/sp.sqrt(Cg.diagonal().prod())
                rho_n[p1,p2] = Cn[0,1]/sp.sqrt(Cn.diagonal().prod())
                rho_g[p2,p1] = rho_g[p1,p2]
                rho_n[p2,p1] = rho_n[p1,p2]
        #3. init
        Cg0 = rho_g*sp.dot(sp.sqrt(var[:,0:1]),sp.sqrt(var[:,0:1].T))
        Cn0 = rho_n*sp.dot(sp.sqrt(var[:,1:2]),sp.sqrt(var[:,1:2].T))
        offset_g = abs(sp.minimum(sp.linalg.eigh(Cg0)[0].min(),0))+1e-4
        offset_n = abs(sp.minimum(sp.linalg.eigh(Cn0)[0].min(),0))+1e-4
        Cg0+=offset_g*sp.eye(self.P)
        Cn0+=offset_n*sp.eye(self.P)
        Lg = sp.linalg.cholesky(Cg0)
        Ln = sp.linalg.cholesky(Cn0)
        Cg_params0 = sp.concatenate([Lg[:,p][:p+1] for p in range(self.P)])
        Cn_params0 = sp.concatenate([Ln[:,p][:p+1] for p in range(self.P)])
        scales0 = sp.concatenate([Cg_params0,1e-2*sp.ones(1),Cn_params0,1e-2*sp.ones(1)])

        return scales0
42143
def plot_best_worst_fits(assignments_df, data, modality_col='Modality',
                         score='$\log_2 K$'):
    """"""
    ncols = 2
    nrows = len(assignments_df.groupby(modality_col).groups.keys())

    fig, axes = plt.subplots(nrows=nrows, ncols=ncols,
                             figsize=(nrows*4, ncols*6))

    axes_iter = axes.flat

    fits = 'Highest', 'Lowest'

    for modality, df in assignments_df.groupby(modality_col):
        df = df.sort_values(score)

        color = MODALITY_TO_COLOR[modality]

        for fit in fits:
            if fit == 'Highest':
                ids = df['Feature ID'][-10:]
            else:
                ids = df['Feature ID'][:10]
            fit_psi = data[ids]
            tidy_fit_psi = fit_psi.stack().reset_index()
            tidy_fit_psi = tidy_fit_psi.rename(columns={'level_0': 'Sample ID',
                                                        'level_1':
                                                            'Feature ID',
                                                        0: '$\Psi$'})
            if tidy_fit_psi.empty:
                continue
            ax = six.next(axes_iter)
            violinplot(x='Feature ID', y='$\Psi$', data=tidy_fit_psi,
                       color=color, ax=ax)
            ax.set(title='{} {} {}'.format(fit, score, modality), xticks=[])
    sns.despine()
    fig.tight_layout()
329466
def convenience_calc_log_likelihood(self, params):
        
        shapes, intercepts, betas = self.convenience_split_params(params)

        args = [betas,
                self.design_3d,
                self.alt_id_vector,
                self.rows_to_obs,
                self.rows_to_alts,
                self.rows_to_mixers,
                self.choice_vector,
                self.utility_transform]

        kwargs = {"ridge": self.ridge, "weights": self.weights}
        log_likelihood = general_log_likelihood(*args, **kwargs)

        return log_likelihood
6915
def _get_fit(self, time: int) -> typing.Tuple[np.ndarray, np.ndarray, float, float]:
        
        rawdata = self.averagedata[:, time]
        domain = np.arange(len(rawdata))
        datalength = len(domain)
        coefficients = np.zeros((datalength, self.function_number + 2))
        coefficients[:, 0] = 1
        coefficients[:, 1] = domain
        for i in range(self.function_number):
            coefficients[:, 2 + i] = self._gaussian_function(datalength, domain, 1, i)
        betas = linalg.inv(coefficients.transpose().dot(coefficients)).dot(coefficients.transpose().dot(rawdata))
        predicted_values = coefficients.dot(betas)
        residuals = rawdata - predicted_values
        error = np.sqrt(residuals.transpose().dot(residuals) / (datalength - (self.function_number + 2)))
        return predicted_values, residuals, residuals.mean(), error
1151312
def sparse_hessian_log_joint(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        # Collect the Gaussian LDS prior terms
        J_diag, J_upper_diag = self.sparse_J_prior
        H_diag, H_upper_diag = -J_diag, -J_upper_diag

        # Collect the likelihood terms
        H_diag += self.hessian_local_log_likelihood(x)

        # Subtract a little bit to ensure negative definiteness
        H_diag -= 1e-8 * np.eye(D)

        return H_diag, H_upper_diag
566131
def log_likelihood(self, x, x_dims=(128, 128), standard=False):
        
        mu, std, h, w = self.params
        shape = tf.stack([32, x_dims[0], x_dims[1], self.num_latent, 1, 1])
        x_reshape = tf.reshape(x, shape)
        c = - 0.5 * math.log(2 * math.pi)
        if standard is False:
            density = c - tf.log(std + 1e-10) - (x_reshape - mu) ** 2 / (2 * std**2)
        else:
            density = c - x_reshape ** 2 / 2
        # sum over all importance weights. average over all eq_samples
        return tf.reduce_mean(tf.reduce_sum(density, axis=2), axis=(1, 2))
1246024
def report_fit(self):
        
        if not self.fitted:
            print('Model not yet fit.')
            return

        print('Null Log-liklihood: {0:.3f}'.format(
            self.log_likelihoods['null']))
        print('Log-liklihood at convergence: {0:.3f}'.format(
            self.log_likelihoods['convergence']))
        print('Log-liklihood Ratio: {0:.3f}\n'.format(
            self.log_likelihoods['ratio']))

        tbl = PrettyTable(
            ['Component', ])
        tbl = PrettyTable()

        tbl.add_column('Component', self.fit_parameters.index.values)
        for col in ('Coefficient', 'Std. Error', 'T-Score'):
            tbl.add_column(col, self.fit_parameters[col].values)

        tbl.align['Component'] = 'l'
        tbl.float_format = '.3'

        print(tbl)
1106937
def convenience_calc_log_likelihood(self, params):
        
        shapes, intercepts, betas = self.convenience_split_params(params)

        args = [betas,
                self.design,
                self.alt_id_vector,
                self.rows_to_obs,
                self.rows_to_alts,
                self.choice_vector,
                self.utility_transform]

        kwargs = {"intercept_params": intercepts,
                  "shape_params": shapes,
                  "ridge": self.ridge,
                  "weights": self.weights}
        log_likelihood = cc.calc_log_likelihood(*args, **kwargs)

        return log_likelihood
707795
def _fit_m(D, a0, logp, tol=1e-7, maxiter=1000):
    ''''''
    N,K = D.shape
    s = a0.sum()

    for i in xrange(maxiter):
        m = a0 / s
        a1 = _ipsi(logp + (m*(psi(a0) - logp)).sum())
        a1 = a1/a1.sum() * s

        if norm(a1 - a0) < tol:
            return a1
        a0 = a1

    raise Exception('Failed to converge after {} iterations, s is {}'
            .format(maxiter, s))
190907
def expected_log_likelihood(self, mus, sigmas, y):
        

        # Flatten the covariance
        T = mus.shape[0]
        D = self.D_in
        sigs_vec = sigmas.reshape((T, D ** 2))

        # Compute the log likelihood of each column
        ll = np.zeros((T, self.D_out))
        for n in range(self.D_out):

            an = self.A[n]

            E_loglmbda = np.dot(mus, an)
            ll[:,n] += y[:,n] * E_loglmbda

            # Vectorized log likelihood calculation
            aa_vec = np.outer(an, an).reshape((D ** 2,))
            ll[:,n] = -np.exp(E_loglmbda + 0.5 * np.dot(sigs_vec, aa_vec))

        return ll