From 1e2f27edc68c064dd5784898518c52d9f9e7e2bc Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Tue, 7 Nov 2017 17:00:59 -0500 Subject: [PATCH 1/3] Fixing numerical issue and adding utility function --- brainiak/eventseg/event.py | 79 +++++++++++++++++++++++++++--------- tests/eventseg/test_event.py | 25 ++++++++++++ 2 files changed, 84 insertions(+), 20 deletions(-) diff --git a/brainiak/eventseg/event.py b/brainiak/eventseg/event.py index 9eb7a9854..e0c12abb5 100644 --- a/brainiak/eventseg/event.py +++ b/brainiak/eventseg/event.py @@ -91,21 +91,6 @@ def __init__(self, n_events=2, self.step_var = step_var self.n_iter = n_iter - # Set up transition matrix, with final sink state - # For transition matrix of this form, the transition probability has - # no impact on the final solution, since all valid paths must take - # the same number of transitions - self.p_start = np.zeros((1, self.n_events + 1)) - self.p_start[0, 0] = 1 - self.P = np.vstack((np.hstack(( - 0.5 * np.diag(np.ones(self.n_events)) - + 0.5 * np.diag(np.ones(self.n_events - 1), 1), - np.append(np.zeros((self.n_events - 1, 1)), [[0.5]], axis=0))), - np.append(np.zeros((1, self.n_events)), [[1]], - axis=1))) - self.p_end = np.zeros((1, self.n_events + 1)) - self.p_end[0, -2] = 1 - def fit(self, X, y=None): """Learn a segmentation on training data @@ -266,31 +251,47 @@ def _forward_backward(self, logprob): log_alpha = np.zeros((t, self.n_events + 1)) log_beta = np.zeros((t, self.n_events + 1)) + # Set up transition matrix, with final sink state + # For transition matrix of this form, the transition probability has + # no impact on the final solution, since all valid paths must take + # the same number of transitions + p_start = np.zeros((1, self.n_events + 1)) + p_start[0, 0] = 1 + p_trans = (self.n_events-1)/t + P = np.vstack((np.hstack(( + (1 - p_trans) * np.diag(np.ones(self.n_events)) + + p_trans * np.diag(np.ones(self.n_events - 1), 1), + np.append(np.zeros((self.n_events - 1, 1)), [[p_trans]], axis=0))), + np.append(np.zeros((1, self.n_events)), [[1]], + axis=1))) + p_end = np.zeros((1, self.n_events + 1)) + p_end[0, -2] = 1 + # Forward pass for t in range(t): if t == 0: - log_alpha[0, :] = self._log(self.p_start) + logprob[0, :] + log_alpha[0, :] = self._log(p_start) + logprob[0, :] else: log_alpha[t, :] = self._log(np.exp(log_alpha[t - 1, :]) - .dot(self.P)) + logprob[t, :] + .dot(P)) + logprob[t, :] log_scale[t] = np.logaddexp.reduce(log_alpha[t, :]) log_alpha[t] -= log_scale[t] # Backward pass - log_beta[-1, :] = self._log(self.p_end) - log_scale[-1] + log_beta[-1, :] = self._log(p_end) - log_scale[-1] for t in reversed(range(t - 1)): obs_weighted = log_beta[t + 1, :] + logprob[t + 1, :] offset = np.max(obs_weighted) log_beta[t, :] = offset + self._log( - np.exp(obs_weighted - offset).dot(self.P.T)) - log_scale[t] + np.exp(obs_weighted - offset).dot(P.T)) - log_scale[t] # Combine and normalize log_gamma = log_alpha + log_beta log_gamma -= np.logaddexp.reduce(log_gamma, axis=1, keepdims=True) ll = np.sum(log_scale[:(t - 1)]) + np.logaddexp.reduce( - log_alpha[-1, :] + log_scale[-1] + self._log(self.p_end), axis=1) + log_alpha[-1, :] + log_scale[-1] + self._log(p_end), axis=1) log_gamma = log_gamma[:, :-1] @@ -401,3 +402,41 @@ def predict(self, X): X = check_array(X) segments, test_ll = self.find_events(X) return np.argmax(segments, axis=1) + + def calc_weighted_event_var(self, D, weights, event_pat): + """Computes normalized weighted variance around event pattern + + Utility function for computing variance in a training set of weighted + event examples. For each event, the sum of squared differences for all + timepoints from the event pattern is computed, and then the weights + specify how much each of these differences contributes to the + variance (normalized by the number of voxels). + + Parameters + ---------- + D : timepoint by voxel ndarray + fMRI data for which to compute event variances + + weights : timepoint by event ndarray + specifies relative weights of timepoints for each event + + event_pat : voxel by event ndarray + mean event patterns to compute variance around + + Returns + ------- + ev_var : ndarray of variances for each event + """ + Dz = stats.zscore(D, axis=1, ddof=1) + ev_var = np.empty(event_pat.shape[1]) + for e in range(event_pat.shape[1]): + # Only compute variances for weights > 0.1% of max weight + nz = weights[:, e] > np.max(weights[:, e])/1000 + sumsq = np.dot(weights[nz, e], + np.sum(np.square(Dz[nz, :] - + event_pat[:, e]), axis=1)) + ev_var[e] = sumsq/(np.sum(weights[nz, e]) - + np.sum(np.square(weights[nz, e])) / + np.sum(weights[nz, e])) + ev_var = ev_var / D.shape[1] + return ev_var diff --git a/tests/eventseg/test_event.py b/tests/eventseg/test_event.py index bfe227d61..fe0f15650 100644 --- a/tests/eventseg/test_event.py +++ b/tests/eventseg/test_event.py @@ -49,3 +49,28 @@ def test_event_transfer(): events = np.argmax(seg, axis=1) assert np.array_equal(events, [0, 0, 0, 1, 1, 1, 1]),\ "Failed to correctly transfer two events to new data" + + +def test_weighted_var(): + es = brainiak.eventseg.event.EventSegment(2) + + D = np.zeros((8, 4)) + for t in range(4): + D[t, :] = (1/np.sqrt(4/3)) * np.array([-1, -1, 1, 1]) + for t in range(4, 8): + D[t, :] = (1 / np.sqrt(4 / 3)) * np.array([1, 1, -1, -1]) + mean_pat = D[[0, 4], :].T + + weights = np.zeros((8, 2)) + weights[:, 0] = [1, 1, 1, 1, 0, 0, 0, 0] + weights[:, 1] = [0, 0, 0, 0, 1, 1, 1, 1] + assert np.array_equal( + es.calc_weighted_event_var(D, weights, mean_pat), [0, 0]),\ + "Failed to compute variance with 0/1 weights" + + weights[:, 0] = [1, 1, 1, 1, 0.5, 0.5, 0.5, 0.5] + weights[:, 1] = [0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1] + true_var = (4 * 0.5 * 12)/(6 - 5/6) * np.ones(2) / 4 + assert np.allclose( + es.calc_weighted_event_var(D, weights, mean_pat), true_var),\ + "Failed to compute variance with fractional weights" From 30b69e889aae62a7506723fb1f583d365f4a34da Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Tue, 7 Nov 2017 17:12:09 -0500 Subject: [PATCH 2/3] Adding news doc --- docs/newsfragments/292.feature | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 docs/newsfragments/292.feature diff --git a/docs/newsfragments/292.feature b/docs/newsfragments/292.feature new file mode 100644 index 000000000..dcbf560af --- /dev/null +++ b/docs/newsfragments/292.feature @@ -0,0 +1,2 @@ +eventseg: Fixed numerical instability bug, added utility function for weighted variance + From 500d124e86eb64fa0327bae22a47c65c8a515dcd Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Tue, 7 Nov 2017 17:22:06 -0500 Subject: [PATCH 3/3] Fixing test to ensure it runs to convergence --- tests/eventseg/test_event.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/eventseg/test_event.py b/tests/eventseg/test_event.py index fe0f15650..be6ef345e 100644 --- a/tests/eventseg/test_event.py +++ b/tests/eventseg/test_event.py @@ -31,7 +31,7 @@ def test_fit_shapes(): def test_simple_boundary(): - es = brainiak.eventseg.event.EventSegment(2, n_iter=10) + es = brainiak.eventseg.event.EventSegment(2) sample_data = np.asarray([[1, 1, 1, 0, 0, 0, 0], [0, 0, 0, 1, 1, 1, 1]]) es.fit(sample_data.T)