diff --git a/__pycache__/pmf_analysis.cpython-37.pyc b/__pycache__/pmf_analysis.cpython-37.pyc
index 5d73cb225dcfed4bac0c68eecdffceec3453a4d1..b3f10b290da8b80978831d35ff6e92f052377bf5 100644
Binary files a/__pycache__/pmf_analysis.cpython-37.pyc and b/__pycache__/pmf_analysis.cpython-37.pyc differ
diff --git a/all_first_pmfs_typeless b/all_first_pmfs_typeless
deleted file mode 100644
index 336970346beed422f9657b5d9db49519023836d9..0000000000000000000000000000000000000000
Binary files a/all_first_pmfs_typeless and /dev/null differ
diff --git a/pmf_analysis.py b/analysis_pmf.py
similarity index 96%
rename from pmf_analysis.py
rename to analysis_pmf.py
index 9e7a0338c09ff15877860996f52d72999c16cd15..b07ba8a8e44f3f518bd19fc3d8e0b2a8aa2073e4 100644
--- a/pmf_analysis.py
+++ b/analysis_pmf.py
@@ -1,6 +1,7 @@
import pickle
import matplotlib.pyplot as plt
import numpy as np
+import seaborn as sns
type2color = {0: 'green', 1: 'blue', 2: 'red'}
@@ -21,7 +22,7 @@ def pmf_type(pmf):
if __name__ == "__main__":
- all_first_pmfs_typeless = pickle.load(open("all_first_pmfs_typeless", 'rb'))
+ all_first_pmfs_typeless = pickle.load(open("all_first_pmfs_typeless.p", 'rb'))
all_pmfs = pickle.load(open("all_pmfs.p", 'rb'))
fewer_states_side = []
@@ -35,6 +36,9 @@ if __name__ == "__main__":
animal_biases[1] += 1
fewer_states_side.append(np.min(animal_biases / animal_biases.sum()))
plt.hist(fewer_states_side)
+ plt.title("Mixed biases")
+ plt.ylabel("# of mice")
+ plt.xlabel("min(% left biased states, % right biased states)")
sns.despine()
plt.tight_layout()
plt.savefig("./meeting_figures/proportion_other_bias")
@@ -79,6 +83,7 @@ if __name__ == "__main__":
plt.axhline(0.5, color='grey')
plt.savefig("./meeting_figures/bias_change_{}".format(total_counter))
plt.close()
+ print("Bias counters")
print(total_counter)
print(bias_counter)
print(tendency_counter)
@@ -96,6 +101,8 @@ if __name__ == "__main__":
# plt.show()
plt.hist(pmf_ranges, bins=40)
plt.title("Type 2 ranges")
+ plt.ylabel("# of type 2 states")
+ plt.xlabel("PMF range")
plt.show()
lapses = []
@@ -105,6 +112,8 @@ if __name__ == "__main__":
lapses.append(max(pmf[0], 1 - pmf[-1]))
plt.hist(lapses, bins=40)
plt.title("Higher lapse rate of type != 0")
+ plt.ylabel("# of non-type-0 states")
+ plt.xlabel("Higher lapse rate")
plt.show()
n_rows, n_cols = 5, 6
@@ -222,8 +231,9 @@ if __name__ == "__main__":
all_changing_pmf_names = pickle.load(open("changing_pmf_names.p", 'rb'))
plt.figure(figsize=(16, 9))
+ plt.title("Type changing PMFs")
for i, (pmf, name) in enumerate(zip(all_changing_pmfs, all_changing_pmf_names)):
- plt.subplot(4, 7, i + 1)
+ plt.subplot(5, 7, i + 1)
plt.title(name)
for p in pmf[1]:
plt.plot(np.where(pmf[0])[0], p[pmf[0]], color=type2color[pmf_type(p)])
@@ -252,7 +262,6 @@ if __name__ == "__main__":
all_first_pmfs = pickle.load(open("pmfs_temp.p", 'rb'))
n_rows, n_cols = 1, 3
_, axs = plt.subplots(n_rows, n_cols, figsize=(16, 9))
- counter = [[0, 0], [0, 0]]
save_title = "all types" if True else "KS014 types"
if save_title == "KS014 types":
all_first_pmfs_typeless = {'KS014': all_first_pmfs_typeless['KS014']}
@@ -279,14 +288,13 @@ if __name__ == "__main__":
axs[0].set_ylabel("P(rightwards)", size=label_size)
axs[0].set_xlabel("Contrasts", size=label_size)
- print(counter)
plt.tight_layout()
plt.savefig(save_title)
plt.show()
if save_title == "KS014 types":
quit()
- # Type 2 PMFs
+ # Type 1 PMF specifications
counter = 0
fig, ax = plt.subplots(1, 3, figsize=(16, 9))
for key in all_first_pmfs_typeless:
diff --git a/analysis_regression.py b/analysis_regression.py
new file mode 100644
index 0000000000000000000000000000000000000000..fe165ae1aded90bb19b0826ec79f7735284d2608
--- /dev/null
+++ b/analysis_regression.py
@@ -0,0 +1,52 @@
+"""Cool function for offsetting points in a scatter."""
+import pickle
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import pearsonr
+import seaborn as sns
+
+fontsize = 22
+
+if __name__ == "__main__":
+ regressions = np.array(pickle.load(open("regressions.p", 'rb')))
+
+ assert (regressions[:, 0] == np.sum(regressions[:, 2:], 1)).all()
+
+ print(pearsonr(regressions[:, 0], regressions[:, 1]))
+
+ offset = 0.125
+ # which x values exist
+ for x in np.unique(regressions[:, 0]):
+ # which and how many ys are associated with this x
+ tempa, tempb = np.unique(regressions[regressions[:, 0] == x, 1], return_counts=True)
+ for y, count in zip(tempa, tempb):
+ # plot them
+ plt.scatter(x + np.linspace(-count + 1, count - 1, count) / 2 * offset, [y] * count, color='grey')
+ plt.ylabel("# of sessions", size=fontsize)
+ plt.xlabel("# of regressed sessions", size=fontsize)
+ sns.despine()
+
+ plt.tight_layout()
+ plt.savefig("Regression vs session length")
+ plt.show()
+
+ # histogram of regressions per type
+ plt.bar([0, 1, 2], np.sum(regressions[:, 2:], 0), color='grey')
+ plt.ylabel("# of regressed sessions", size=fontsize)
+ plt.gca().set_xticks([0, 1, 2], ['Type 1', 'Type 2', 'Type 3'], size=fontsize)
+
+ sns.despine()
+ plt.tight_layout()
+ plt.savefig("# of regressions per type")
+ plt.show()
+
+ # histogram of number of mice with regressions in the different types
+ num_of_mice = [*np.sum(regressions[:, 2:] >= 1, 0), np.sum(regressions[:, 0] >= 1)]
+ plt.bar([0, 1, 2, 3], num_of_mice, color='grey')
+ plt.ylabel("# of mice with regressions", size=fontsize)
+ plt.gca().set_xticks([0, 1, 2, 3], ['Type 1', 'Type 2', 'Type 3', 'Any'], size=fontsize)
+
+ sns.despine()
+ plt.tight_layout()
+ plt.savefig("# of mice with regressions per type")
+ plt.show()
diff --git a/analysis_state_intros.py b/analysis_state_intros.py
new file mode 100644
index 0000000000000000000000000000000000000000..c00985af68fa992ce0fba0a42df16274370b7fa9
--- /dev/null
+++ b/analysis_state_intros.py
@@ -0,0 +1,55 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import pickle
+
+fontsize = 16
+
+def type_hist(data, title=''):
+ highest = int(data.max())
+ if (data % 1 == 0).all():
+ bins = np.arange(highest + 2) - 0.5
+ else:
+ bins = np.histogram(data)[1]
+ hist_max = 0
+ for i in range(data.shape[1]):
+ hist_max = max(hist_max, np.histogram(data[:, i], bins)[0].max())
+
+ plt.subplot(3, 1, 1)
+
+ if title != '':
+ plt.title(title, size=fontsize + 4)
+
+ assert np.histogram(data[:, 0])[0].sum() == np.histogram(data[:, 0], bins)[0].sum()
+ plt.hist(data[:, 0], alpha=1/3, label="type 1", align='mid', bins=bins, color='grey')
+ plt.xlim(-0.5, highest + 1)
+ plt.ylim(0, hist_max + 1)
+ plt.ylabel("Type 1", size=fontsize)
+
+ plt.subplot(3, 1, 2)
+ assert np.histogram(data[:, 1])[0].sum() == np.histogram(data[:, 1], bins)[0].sum()
+ plt.hist(data[:, 1], alpha=1/3, label="type 2", align='mid', bins=bins, color='grey')
+ plt.xlim(-0.5, highest + 1)
+ plt.ylim(0, hist_max + 1)
+ plt.ylabel("Type 2", size=fontsize)
+
+ plt.subplot(3, 1, 3)
+ assert np.histogram(data[:, 2])[0].sum() == np.histogram(data[:, 2], bins)[0].sum()
+ plt.hist(data[:, 2], alpha=1/3, label="type 3", align='mid', bins=bins, color='grey')
+ plt.xlim(-0.5, highest + 1)
+ plt.ylim(0, hist_max + 1)
+ plt.ylabel("Type 3", size=fontsize)
+
+ # plt.legend()
+ plt.show()
+
+
+if __name__ == "__main__":
+ all_intros = pickle.load(open("all_intros.p", 'rb'))
+ all_intros_div = pickle.load(open("all_intros_div.p", 'rb'))
+ all_states_per_type = pickle.load(open("all_states_per_type.p", 'rb'))
+
+ # There are 5 mice with 0 type 2 intros, but only 3 mice with no type 2 stats.
+ # That is because they come up, but don't explain the necessary 50% to start phase 2.
+ type_hist(all_intros, title='State introductions')
+ type_hist(all_intros_div, title='Session-normalised state introductions')
+ type_hist(all_states_per_type, title='Average states used per session')
diff --git a/dyn_glm_chain_analysis.py b/dyn_glm_chain_analysis.py
index c44c7897447e7df141ad2fd4e516f65052f99d28..4924143b4f44716c903d7e798a60125d432c262e 100644
--- a/dyn_glm_chain_analysis.py
+++ b/dyn_glm_chain_analysis.py
@@ -606,7 +606,7 @@ def pmf_regressions(states_by_session, pmfs, durs):
if current_best_state == -1 or state_perfs[current_best_state] < state_perfs[s]:
current_best_state = s
- if state_perfs[np.argmax(states_by_session[:, sess])] + 0.05 < state_perfs[current_best_state] and state_counter[np.argmax(states_by_session[:, sess])] > 1:
+ if state_perfs[np.argmax(states_by_session[:, sess])] + 0.025 < state_perfs[current_best_state] and state_counter[np.argmax(states_by_session[:, sess])] > 1:
counter += 1
if sess < durs[0]:
types[0] += 1
@@ -1429,51 +1429,6 @@ def compare_performance(cnas, contrasts=(1, 1), title=""):
plt.show()
-def type_hist(data):
- highest = int(data.max())
- if (data % 1 == 0).all():
- bins = np.arange(highest + 2) - 0.5
- else:
- bins = np.histogram(data)[1]
- hist_max = 0
- for i in range(data.shape[1]):
- hist_max = max(hist_max, np.histogram(data[:, i], bins)[0].max())
-
- plt.subplot(3, 1, 1)
- assert np.histogram(data[:, 0])[0].sum() == np.histogram(data[:, 0], bins)[0].sum()
- plt.hist(data[:, 0], alpha=1/3, label="type 1", align='mid', bins=bins)
- plt.xlim(-0.5, highest + 1)
- plt.ylim(0, hist_max + 1)
-
- plt.subplot(3, 1, 2)
- assert np.histogram(data[:, 1])[0].sum() == np.histogram(data[:, 1], bins)[0].sum()
- plt.hist(data[:, 1], alpha=1/3, label="type 2", align='mid', bins=bins)
- plt.xlim(-0.5, highest + 1)
- plt.ylim(0, hist_max + 1)
-
- plt.subplot(3, 1, 3)
- assert np.histogram(data[:, 2])[0].sum() == np.histogram(data[:, 2], bins)[0].sum()
- plt.hist(data[:, 2], alpha=1/3, label="type 3", align='mid', bins=bins)
- plt.xlim(-0.5, highest + 1)
- plt.ylim(0, hist_max + 1)
-
- # plt.legend()
- plt.show()
-
-
-if 0:
- all_intros = pickle.load(open("all_intros.p", 'rb'))
- all_intros_div = pickle.load(open("all_intros_div.p", 'rb'))
- all_states_per_type = pickle.load(open("all_states_per_type.p", 'rb'))
-
- # There are 5 mice with 0 type 2 intros, but only 3 mice with no type 2 stats.
- # That is because they come up, but don't explain the necessary 50% to start phase 2.
- type_hist(all_intros)
- type_hist(all_intros_div)
- type_hist(all_states_per_type)
- quit()
-
-
if __name__ == "__main__":
fit_type = ['prebias', 'bias', 'all', 'prebias_plus', 'zoe_style'][0]
@@ -1485,7 +1440,7 @@ if __name__ == "__main__":
r_hats = json.load(open("canonical_info_r_hats.json", 'r'))
no_good_pcas = ['NYU-06', 'SWC_023']
subjects = list(loading_info.keys())
- # subjects = ['ZM_1897']
+ subjects = ['KS014']
# meh pmfs: KS021
print(subjects)
@@ -1533,7 +1488,6 @@ if __name__ == "__main__":
print(subject)
try:
- continue
test = pickle.load(open("multi_chain_saves/canonical_result_{}_{}.p".format(subject, fit_type), 'rb'))
print('loaded canoncial result')
@@ -1548,7 +1502,7 @@ if __name__ == "__main__":
states, pmfs, durs, _, contrast_intro_type, intros_by_type, undiv_intros, states_per_type = state_development(test, [s for s in state_sets if len(s) > 40], mode_indices, show=0, separate_pmf=1, type_coloring=True)
regression = pmf_regressions(states, pmfs, durs)
regressions.append(regression)
- continue
+ quit()
# compare_pmfs(test, [s for s in state_sets if len(s) > 40], mode_indices, [4, 5], states, pmfs, title="{} convergence pmf".format(subject))
# compare_weights(test, [s for s in state_sets if len(s) > 40], mode_indices, [4, 5], states, title="{} convergence weights".format(subject))
# quit()
@@ -1774,6 +1728,7 @@ if __name__ == "__main__":
continue
else:
print("Making canonical result")
+ results = []
for j, (seed, fit_num) in enumerate(zip(loading_info[subject]['seeds'], loading_info[subject]['fit_nums'])):
if j in loading_info[subject]['ignore']:
continue
@@ -1849,8 +1804,8 @@ if __name__ == "__main__":
print("Correlations")
from scipy.stats import pearsonr
print(pearsonr(abs_state_durs[:, 0], abs_state_durs[:, 1]))
- print(pearsonr(abs_state_durs[:, 2], abs_state_durs[:, 1]))
print(pearsonr(abs_state_durs[:, 0], abs_state_durs[:, 2]))
+ print(pearsonr(abs_state_durs[:, 1], abs_state_durs[:, 2]))
print(pearsonr(abs_state_durs.sum(1), abs_state_durs[:, 0]))
# (0.7338297529946006, 2.6332570579118393e-06)