diff --git a/src/rsatoolbox/data/noise.py b/src/rsatoolbox/data/noise.py index ebc90eac..3220afc1 100755 --- a/src/rsatoolbox/data/noise.py +++ b/src/rsatoolbox/data/noise.py @@ -32,16 +32,18 @@ def _check_demean(matrix): matrix = matrix - np.mean(matrix, axis=0, keepdims=True) dof = matrix.shape[0] - 1 elif matrix.ndim == 3: - matrix -= np.mean(matrix, axis=2, keepdims=True) - dof = (matrix.shape[0] - 1) * matrix.shape[2] + matrix -= np.nanmean(matrix, axis=2, keepdims=True) + # dof = (matrix.shape[0] - 1) * matrix.shape[2] JohnMark checking + dof = matrix.shape[0] * (matrix.shape[2] - 1) matrix = matrix.transpose(0, 2, 1).reshape( matrix.shape[0] * matrix.shape[2], matrix.shape[1]) + matrix = matrix[~np.isnan(matrix).any(axis=1)] else: raise ValueError('Matrix for covariance estimation has wrong # of dimensions!') return matrix, dof -def _estimate_covariance(matrix, dof, method): +def _estimate_covariance(matrix, dof, method, lamb_opt=None): """ calls the right covariance estimation function based on the ""method" argument Args: @@ -64,9 +66,9 @@ def _estimate_covariance(matrix, dof, method): dof = dof_nat # calculate sample covariance matrix s if method == 'shrinkage_eye': - cov_mat = _covariance_eye(matrix, dof) + cov_mat = _covariance_eye(matrix, dof, lamb_opt) elif method == 'shrinkage_diag': - cov_mat = _covariance_diag(matrix, dof) + cov_mat = _covariance_diag(matrix, dof, lamb_opt) elif method == 'diag': cov_mat = _variance(matrix, dof) elif method == 'full': @@ -108,7 +110,7 @@ def _covariance_full(matrix, dof): return np.einsum('ij, ik-> jk', matrix, matrix, optimize=True) / dof -def _covariance_eye(matrix, dof): +def _covariance_eye(matrix, dof, lamb_opt): """ computes the sample covariance matrix from a 2d-array. matrix should be demeaned before! @@ -143,17 +145,21 @@ def _covariance_eye(matrix, dof): # calculate the scalar estimators to find the optimal shrinkage: # m, d^2, b^2 as in Ledoit & Wolfe paper m = np.sum(np.diag(s)) / s.shape[0] - d2 = np.sum((s - m * np.eye(s.shape[0])) ** 2) - b2 = min(d2, b2) - # shrink covariance matrix - s_shrink = b2 / d2 * m * np.eye(s.shape[0]) \ - + (d2-b2) / d2 * s + if lamb_opt is None: + d2 = np.sum((s - m * np.eye(s.shape[0])) ** 2) + b2 = min(d2, b2) + # shrink covariance matrix + s_shrink = b2 / d2 * m * np.eye(s.shape[0]) \ + + (d2 - b2) / d2 * s + else: + s_shrink = lamb_opt * m * np.eye(s.shape[0]) \ + + (1 - lamb_opt) * s # correction for degrees of freedom s_shrink = s_shrink * matrix.shape[0] / dof return s_shrink -def _covariance_diag(matrix, dof, mem_threshold=(10**9)/8): +def _covariance_diag(matrix, dof, lamb_opt=None, mem_threshold=(10 ** 9) / 8): """ computes the sample covariance matrix from a 2d-array. matrix should be demeaned before! @@ -186,15 +192,29 @@ def _covariance_diag(matrix, dof, mem_threshold=(10**9)/8): s = s_sum / dof var = np.diag(s) std = np.sqrt(var) - s_mean = s_sum / np.expand_dims(std, 0) / np.expand_dims(std, 1) / (matrix.shape[0] - 1) - s2_mean = s2_sum / np.expand_dims(var, 0) / np.expand_dims(var, 1) / (matrix.shape[0] - 1) - var_hat = matrix.shape[0] / dof ** 2 \ - * (s2_mean - s_mean ** 2) mask = ~np.eye(s.shape[0], dtype=bool) - lamb = np.sum(var_hat[mask]) / np.sum(s_mean[mask] ** 2) - lamb = max(min(lamb, 1), 0) - scaling = np.eye(s.shape[0]) + (1-lamb) * mask + if lamb_opt is None: + # s_mean = s_sum / np.expand_dims(std, 0) / np.expand_dims(std, 1) / (matrix.shape[0] - 1) JohnMark check + # s2_mean = s2_sum / np.expand_dims(var, 0) / np.expand_dims(var, 1) / (matrix.shape[0] - 1) + s_mean = s_sum / np.expand_dims(std, 0) / np.expand_dims(std, 1) / dof + s2_mean = s2_sum / np.expand_dims(var, 0) / np.expand_dims(var, 1) / dof + var_hat = matrix.shape[0] / dof ** 2 \ + * (s2_mean - s_mean ** 2) + lamb = np.sum(var_hat[mask]) / np.sum(s_mean[mask] ** 2) + lamb = max(min(lamb, 1), 0) + else: + lamb = lamb_opt + scaling = np.eye(s.shape[0]) + (1 - lamb) * mask s_shrink = s * scaling + + # mean_shrunk_eigenvalue = np.mean(np.linalg.eigvals(s_shrink)) + # mean_full_eigenvalue = np.mean(np.linalg.eigvals(s)) + # print(f"data shape: {matrix.shape}") + # print(f"dof: {dof}") + # print(f"lambda: {lamb}") + # print(f"mean var: {np.mean(var)}") + # print(f"mean full eigenvalue: {mean_full_eigenvalue}") + # print(f"mean shrunk eigenvalue: {mean_shrunk_eigenvalue}") return s_shrink @@ -272,7 +292,7 @@ def prec_from_residuals(residuals, dof=None, method='shrinkage_diag'): return prec -def cov_from_measurements(dataset, obs_desc, dof=None, method='shrinkage_diag'): +def cov_from_measurements(dataset, obs_desc, dof=None, method='shrinkage_diag', lamb_opt=None): """ Estimates a covariance matrix from measurements. Allows for shrinkage estimates. Use 'method' to choose which estimation method is used. @@ -311,11 +331,11 @@ def cov_from_measurements(dataset, obs_desc, dof=None, method='shrinkage_diag'): "obs_desc not contained in the dataset's obs_descriptors" tensor, _ = dataset.get_measurements_tensor(obs_desc) # calculate sample covariance matrix s - cov_mat = _estimate_covariance(tensor, dof, method) + cov_mat = _estimate_covariance(tensor, dof, method, lamb_opt) return cov_mat -def prec_from_measurements(dataset, obs_desc, dof=None, method='shrinkage_diag'): +def prec_from_measurements(dataset, obs_desc, dof=None, method='shrinkage_diag', lamb_opt=None): """ Estimates the covariance matrix from measurements and finds its multiplicative inverse (= the precision matrix) @@ -337,7 +357,7 @@ def prec_from_measurements(dataset, obs_desc, dof=None, method='shrinkage_diag') numpy.ndarray (or list): sigma_p: precision matrix over channels """ - cov = cov_from_measurements(dataset, obs_desc, dof=dof, method=method) + cov = cov_from_measurements(dataset, obs_desc, dof=dof, method=method, lamb_opt=lamb_opt) if not isinstance(cov, np.ndarray): prec = [None] * len(cov) for i, cov_i in enumerate(cov): diff --git a/src/rsatoolbox/rdm/calc.py b/src/rsatoolbox/rdm/calc.py index 227c1635..1e607c48 100644 --- a/src/rsatoolbox/rdm/calc.py +++ b/src/rsatoolbox/rdm/calc.py @@ -94,6 +94,12 @@ def calc_rdm( elif method == 'crossnobis': rdm = calc_rdm_crossnobis(dataset, descriptor, noise, cv_descriptor, remove_mean) + elif method == 'dotproduct': + rdm = calc_rdm_dotproduct(dataset, descriptor) + elif method == 'mean_profile': + rdm = calc_rdm_mean_profile(dataset, descriptor) + elif method == 'norm_profile': + rdm = calc_rdm_norm_profile(dataset, descriptor) elif method == 'poisson': rdm = calc_rdm_poisson(dataset, descriptor, prior_lambda=prior_lambda, @@ -105,8 +111,12 @@ def calc_rdm( prior_weight=prior_weight) else: raise NotImplementedError - if descriptor is not None: + if (descriptor is not None) and (method not in ['mean_profile', 'norm_profile']): rdm.sort_by(**{descriptor: 'alpha'}) + else: + desc = np.unique(np.array(dataset.obs_descriptors[descriptor])) + inds = desc.argsort() + rdm = rdm[inds] return rdm @@ -209,9 +219,9 @@ def calc_rdm_euclidean( rsatoolbox.rdm.rdms.RDMs: RDMs object with the one RDM """ measurements, desc = _parse_input(dataset, descriptor, remove_mean) - sum_sq_measurements = np.sum(measurements**2, axis=1, keepdims=True) + sum_sq_measurements = np.sum(measurements ** 2, axis=1, keepdims=True) rdm = sum_sq_measurements + sum_sq_measurements.T \ - - 2 * np.dot(measurements, measurements.T) + - 2 * np.dot(measurements, measurements.T) rdm = _extract_triu_(rdm) / measurements.shape[1] return _build_rdms(rdm, dataset, 'squared euclidean', descriptor, desc) @@ -269,7 +279,7 @@ def calc_rdm_mahalanobis(dataset, descriptor=None, noise=None, remove_mean: bool noise = _check_noise(noise, dataset.n_channel) kernel = measurements @ noise @ measurements.T rdm = np.expand_dims(np.diag(kernel), 0) + \ - np.expand_dims(np.diag(kernel), 1) - 2 * kernel + np.expand_dims(np.diag(kernel), 1) - 2 * kernel rdm = _extract_triu_(rdm) / measurements.shape[1] return _build_rdms( rdm, @@ -281,6 +291,71 @@ def calc_rdm_mahalanobis(dataset, descriptor=None, noise=None, remove_mean: bool ) +def calc_rdm_dotproduct( + dataset: DatasetBase, + descriptor: Optional[str] = None, + remove_mean: bool = False): + """ + Args: + dataset (rsatoolbox.data.DatasetBase): + The dataset the RDM is computed from + descriptor (String): + obs_descriptor used to define the rows/columns of the RDM + defaults to one row/column per row in the dataset + remove_mean (bool): + whether the mean of each pattern shall be removed + before calculating dotproducts. + Returns: + rsatoolbox.rdm.rdms.RDMs: RDMs object with the one RDM + """ + measurements, desc = _parse_input(dataset, descriptor, remove_mean) + rdm = measurements @ measurements.T + rdm = _extract_triu_(rdm) + return _build_rdms(rdm, dataset, 'dotproduct', descriptor, desc) + + +def calc_rdm_mean_profile( + dataset: DatasetBase, + descriptor: Optional[str] = None): + """ + Args: + dataset (rsatoolbox.data.DatasetBase): + The dataset the RDM is computed from + descriptor (String): + obs_descriptor used to define the rows/columns of the RDM + defaults to one row/column per row in the dataset + remove_mean (bool): + whether the mean of each pattern shall be removed + before calculating dotproducts. + Returns: + rsatoolbox.rdm.rdms.RDMs: RDMs object with the one RDM + """ + measurements, desc = _parse_input(dataset, descriptor, remove_mean=False) + measurements = measurements.mean(axis=1) + return measurements + + +def calc_rdm_norm_profile( + dataset: DatasetBase, + descriptor: Optional[str] = None): + """ + Args: + dataset (rsatoolbox.data.DatasetBase): + The dataset the RDM is computed from + descriptor (String): + obs_descriptor used to define the rows/columns of the RDM + defaults to one row/column per row in the dataset + remove_mean (bool): + whether the mean of each pattern shall be removed + before calculating dotproducts. + Returns: + rsatoolbox.rdm.rdms.RDMs: RDMs object with the one RDM + """ + measurements, desc = _parse_input(dataset, descriptor, remove_mean=False) + measurements = np.linalg.norm(measurements, axis=1) + return measurements + + def calc_rdm_crossnobis(dataset, descriptor, noise=None, cv_descriptor=None, remove_mean: bool = False): """ @@ -371,7 +446,7 @@ def calc_rdm_crossnobis(dataset, descriptor, noise=None, measurements[i_fold], measurements[j_fold], np.linalg.inv( (variances[i_fold] + variances[j_fold]) / 2) - ) + ) rdms.append(rdm) rdms = np.array(rdms) rdm = np.einsum('ij->j', rdms) / rdms.shape[0] @@ -406,10 +481,10 @@ def calc_rdm_poisson(dataset, descriptor=None, prior_lambda=1, """ measurements, desc = _parse_input(dataset, descriptor) measurements = (measurements + prior_lambda * prior_weight) \ - / (1 + prior_weight) + / (1 + prior_weight) kernel = measurements @ np.log(measurements).T rdm = np.expand_dims(np.diag(kernel), 0) + \ - np.expand_dims(np.diag(kernel), 1) - kernel - kernel.T + np.expand_dims(np.diag(kernel), 1) - kernel - kernel.T rdm = _extract_triu_(rdm) / measurements.shape[1] return _build_rdms(rdm, dataset, 'poisson', descriptor, desc) @@ -455,13 +530,13 @@ def calc_rdm_poisson_cv(dataset, descriptor=None, prior_lambda=1, measurements_test, _, _ = average_dataset_by(data_test, descriptor) measurements_train = (measurements_train + prior_lambda * prior_weight) \ - / (1 + prior_weight) + / (1 + prior_weight) measurements_test = (measurements_test + prior_lambda * prior_weight) \ - / (1 + prior_weight) + / (1 + prior_weight) kernel = measurements_train @ np.log(measurements_test).T rdm = np.expand_dims(np.diag(kernel), 0) + \ - np.expand_dims(np.diag(kernel), 1) - kernel - kernel.T + np.expand_dims(np.diag(kernel), 1) - kernel - kernel.T rdm = _extract_triu_(rdm) / measurements_train.shape[1] return _build_rdms(rdm, dataset, 'poisson_cv', descriptor) @@ -469,7 +544,7 @@ def calc_rdm_poisson_cv(dataset, descriptor=None, prior_lambda=1, def _calc_rdm_crossnobis_single(meas1, meas2, noise) -> NDArray: kernel = meas1 @ noise @ meas2.T rdm = np.expand_dims(np.diag(kernel), 0) + \ - np.expand_dims(np.diag(kernel), 1) - kernel - kernel.T + np.expand_dims(np.diag(kernel), 1) - kernel - kernel.T return _extract_triu_(rdm) / meas1.shape[1] @@ -481,8 +556,8 @@ def _gen_default_cv_descriptor(dataset, descriptor) -> np.ndarray: desc = dataset.obs_descriptors[descriptor] values, counts = np.unique(desc, return_counts=True) assert np.all(counts == counts[0]), ( - 'cv_descriptor generation failed:\n' - + 'different number of observations per pattern') + 'cv_descriptor generation failed:\n' + + 'different number of observations per pattern') n_repeats = counts[0] cv_descriptor = np.zeros_like(desc) for i_val in values: @@ -491,10 +566,10 @@ def _gen_default_cv_descriptor(dataset, descriptor) -> np.ndarray: def _parse_input( - dataset: DatasetBase, - descriptor: Optional[str], - remove_mean: bool = False - ) -> Tuple[np.ndarray, Optional[np.ndarray]]: + dataset: DatasetBase, + descriptor: Optional[str], + remove_mean: bool = False +) -> Tuple[np.ndarray, Optional[np.ndarray]]: if descriptor is None: measurements = dataset.measurements desc = None diff --git a/src/rsatoolbox/rdm/compare.py b/src/rsatoolbox/rdm/compare.py index 9a3e1d55..dea95024 100644 --- a/src/rsatoolbox/rdm/compare.py +++ b/src/rsatoolbox/rdm/compare.py @@ -3,12 +3,16 @@ """ Comparison methods for comparing two RDMs objects """ +import itertools as it import numpy as np +import os +import pickle import scipy.stats from scipy import linalg from scipy.optimize import minimize from scipy.stats._stats import _kendall_dis from scipy.spatial.distance import squareform +from statistics import mode from rsatoolbox.util.matrix import pairwise_contrast_sparse from rsatoolbox.util.matrix import pairwise_contrast from rsatoolbox.util.rdm_utils import _get_n_from_reduced_vectors @@ -132,7 +136,11 @@ def compare_cosine_cov_weighted(rdm1, rdm2, sigma_k=None): """ vector1, vector2, nan_idx = _parse_input_rdms(rdm1, rdm2) - sim = _cosine_cov_weighted(vector1, vector2, sigma_k, nan_idx) + stim_nums = np.array(rdm1.pattern_descriptors['stimulus']) + frozen_inds = list(np.where(stim_nums == -1)[0]) + list(np.where(stim_nums == -2)[0]) + # sigma_k[frozen_inds, :] = 0 + # sigma_k[:, frozen_inds] = 0 + sim = _cosine_cov_weighted(vector1, vector2, frozen_inds, sigma_k, nan_idx) return sim @@ -152,10 +160,14 @@ def compare_correlation_cov_weighted(rdm1, rdm2, sigma_k=None): """ vector1, vector2, nan_idx = _parse_input_rdms(rdm1, rdm2) + stim_nums = np.array(rdm1.pattern_descriptors['stimulus']) + frozen_inds = list(np.where(stim_nums == -1)[0]) + list(np.where(stim_nums == -2)[0]) + # sigma_k[frozen_inds, :] = 0 + # sigma_k[:, frozen_inds] = 0 # compute by subtracting the mean and then calculating cosine similarity vector1 = vector1 - np.mean(vector1, 1, keepdims=True) vector2 = vector2 - np.mean(vector2, 1, keepdims=True) - sim = _cosine_cov_weighted(vector1, vector2, sigma_k, nan_idx) + sim = _cosine_cov_weighted(vector1, vector2, frozen_inds, sigma_k, nan_idx) return sim @@ -259,16 +271,16 @@ def compare_neg_riemannian_distance(rdm1, rdm2, sigma_k=None): n_cond = _get_n_from_length(vector1.shape[1]) if sigma_k is None: sigma_k = np.eye(n_cond) - P = np.block([-1*np.ones((n_cond - 1, 1)), np.eye(n_cond - 1)]) - sigma_k_hat = P@sigma_k@P.T + P = np.block([-1 * np.ones((n_cond - 1, 1)), np.eye(n_cond - 1)]) + sigma_k_hat = P @ sigma_k @ P.T # construct RDM to 2nd-moment (G) transformation - pairs = pairwise_contrast(np.arange(n_cond-1)) + pairs = pairwise_contrast(np.arange(n_cond - 1)) pairs[pairs == -1] = 1 T = np.block([ - [np.eye(n_cond - 1), np.zeros((n_cond-1, vector1.shape[1] - n_cond + 1))], + [np.eye(n_cond - 1), np.zeros((n_cond - 1, vector1.shape[1] - n_cond + 1))], [0.5 * pairs, np.diag(-0.5 * np.ones(vector1.shape[1] - n_cond + 1))]]) - vec_G1 = vector1@np.transpose(T) - vec_G2 = vector2@np.transpose(T) + vec_G1 = vector1 @ np.transpose(T) + vec_G2 = vector2 @ np.transpose(T) sim = _all_combinations(vec_G1, vec_G2, _riemannian_distance, sigma_k_hat) return sim @@ -301,7 +313,7 @@ def _all_combinations(vectors1, vectors2, func, *args, **kwargs): return value -def _cosine_cov_weighted_slow(vector1, vector2, sigma_k=None, nan_idx=None): +def _cosine_cov_weighted_slow(vector1, vector2, frozen_inds=[], sigma_k=None, nan_idx=None): """computes the cosine similarities between two sets of vectors after whitening by their covariance. @@ -327,6 +339,14 @@ def _cosine_cov_weighted_slow(vector1, vector2, sigma_k=None, nan_idx=None): else: n_cond = _get_n_from_reduced_vectors(vector1) v = _get_v(n_cond, sigma_k) + # Now adjust v to account for any frozen patterns. + # v = _correct_covariance_for_frozen_patterns(v, n_cond, frozen_inds) + # Omit any all-zero rows and columns, keeping as a sparse matrix. + nonzero_rows = np.where(v.sum(axis=1) != 0)[0] + v = v[nonzero_rows][:, nonzero_rows] + vector1 = vector1[:, nonzero_rows] + vector2 = vector2[:, nonzero_rows] + # compute V^-1 vector1/2 for all vectors by solving Vx = vector1/2 vector1_m = np.array([scipy.sparse.linalg.cg(v, vector1[i], atol=0)[0] for i in range(vector1.shape[0])]) @@ -343,7 +363,46 @@ def _cosine_cov_weighted_slow(vector1, vector2, sigma_k=None, nan_idx=None): return cos -def _cosine_cov_weighted(vector1, vector2, sigma_k=None, nan_idx=None): +def _correct_covariance_for_frozen_patterns(v, n_cond, frozen_inds): + '''Rules: + 1. If a shared pattern is frozen, set covariance to zero + 2. If one pattern in a pair is frozen, set variance to half. + 3. If both patterns frozen, set variance to zero. + ''' + frozen_inds = set(frozen_inds) + if len(frozen_inds) == 0: + return v + + # Fetch cached value if it exists. + fname = "_".join(str(num) for num in sorted(frozen_inds)) + f"_{n_cond}conds" + "_cov_matrix.pkl" + if os.path.exists(fname): + with open(fname, 'rb') as f: + v = pickle.load(f) + return v + + for (i, (x1, y1)), (j, (x2, y2)) in it.product(enumerate(it.combinations(range(n_cond), 2)), + enumerate(it.combinations(range(n_cond), 2))): + if len(np.unique([x1, y1, x2, y2])) == 4: # if no shared patterns, skip + continue + + if len(np.unique([x1, y1, x2, y2])) == 2: # if the same pair: + num_frozen = np.sum([num in frozen_inds for num in [x1, y1]]) + if num_frozen == 2: + v[i, j] = 0 + elif num_frozen == 1: + v[i, j] /= 2 + + if len(np.unique([x1, y1, x2, y2])) == 3: # if one shared pattern + # Check if the shared pattern is frozen, and if so, set the entry to zero. + rep_ind = mode([x1, y1, x2, y2]) + if rep_ind in frozen_inds: + v[i, j] = 0 + with open(fname, 'wb') as f: + pickle.dump(v, f, pickle.HIGHEST_PROTOCOL) + return v + + +def _cosine_cov_weighted(vector1, vector2, frozen_inds=[], sigma_k=None, nan_idx=None): """computes the cosine angles between two sets of vectors weighted by the covariance If no covariance is given this is computed using the linear CKA, @@ -365,7 +424,7 @@ def _cosine_cov_weighted(vector1, vector2, sigma_k=None, nan_idx=None): """ if (sigma_k is not None) and (sigma_k.ndim >= 2): cos = _cosine_cov_weighted_slow( - vector1, vector2, sigma_k=sigma_k, nan_idx=nan_idx) + vector1, vector2, frozen_inds=frozen_inds, sigma_k=sigma_k, nan_idx=nan_idx) else: if nan_idx is None: nan_idx = np.ones(vector1[0].shape, bool) @@ -428,8 +487,8 @@ def _cov_weighting(vector, nan_idx, sigma_k=None): diag = np.concatenate((np.ones((n_dist, 1)) / 2, np.ones((n_cond, 1)))) # one line version much faster here! vector_w = vector_w - ( - vector_w - @ sumI @ np.linalg.inv(sumI.T @ (diag * sumI)) @ (diag * sumI).T) + vector_w + @ sumI @ np.linalg.inv(sumI.T @ (diag * sumI)) @ (diag * sumI).T) if sigma_k is not None: if sigma_k.ndim == 1: sigma_k_sqrt = np.sqrt(sigma_k) @@ -492,12 +551,13 @@ def _riemannian_distance(vec_G1, vec_G2, sigma_k): negative riemannian distance """ n_cond = _get_n_from_length(len(vec_G1)) - G1 = np.diag(vec_G1[0:(n_cond-1)])+squareform(vec_G1[(n_cond-1):len(vec_G1)]) - G2 = np.diag(vec_G2[0:(n_cond-1)])+squareform(vec_G2[(n_cond-1):len(vec_G2)]) + G1 = np.diag(vec_G1[0:(n_cond - 1)]) + squareform(vec_G1[(n_cond - 1):len(vec_G1)]) + G2 = np.diag(vec_G2[0:(n_cond - 1)]) + squareform(vec_G2[(n_cond - 1):len(vec_G2)]) def fun(theta): return np.sqrt((np.log(linalg.eigvalsh( - np.exp(theta[0]) * G1 + np.exp(theta[1]) * sigma_k, G2))**2).sum()) + np.exp(theta[0]) * G1 + np.exp(theta[1]) * sigma_k, G2)) ** 2).sum()) + theta = minimize(fun, (0, 0), method='Nelder-Mead') neg_riem = -1 * theta.fun return neg_riem @@ -542,8 +602,8 @@ def _tau_a(vector1, vector2): (vector2[1:] != vector2[:-1]), True] cnt = np.diff(np.nonzero(obs)[0]).astype('int64', copy=False) ntie = (cnt * (cnt - 1) // 2).sum() # joint ties - xtie, _, _ = _count_rank_tie(vector1) # ties in x, stats - ytie, _, _ = _count_rank_tie(vector2) # ties in y, stats + xtie, _, _ = _count_rank_tie(vector1) # ties in x, stats + ytie, _, _ = _count_rank_tie(vector2) # ties in y, stats tot = (size * (size - 1)) // 2 # Note that tot = con + dis + (xtie - ntie) + (ytie - ntie) + ntie # = con + dis + xtie + ytie - ntie @@ -569,7 +629,7 @@ def _count_rank_tie(ranks): cnt = cnt[cnt > 1] return ((cnt * (cnt - 1) // 2).sum(), (cnt * (cnt - 1.) * (cnt - 2)).sum(), - (cnt * (cnt - 1.) * (2*cnt + 5)).sum()) + (cnt * (cnt - 1.) * (2 * cnt + 5)).sum()) def _get_v(n_cond, sigma_k):