In [39]:
import numpy as np
from scipy.io import loadmat
from os.path import join
import glob
import pandas as pd
import tensorflow as tf
from sklearn import svm
from sklearn.metrics import accuracy_score
In [59]:
def center_normalize(x):
x_mean_rows = np.mean(x,1).reshape(x.shape[0],1)
x_std_rows = np.std(x,1).reshape(x.shape[0],1)
return (x - x_mean_rows) / x_std_rows
def onehot(y):
ynp=np.array(y)
y_onehot=[0]*len(ynp)
for i,j in enumerate(ynp):
y_onehot[i]=[0]*ynp.max()
y_onehot[i][j-1]=1
return y_onehot
def model_svm(x_train,y_train,x_test,y_test):
C = 1 # SVM regularization parameter
svc = svm.LinearSVC(C=C).fit(x_train, y_train)
y_pred = svc.predict(x_test)
return accuracy_score(y_test, y_pred)
def model_tf_regression(x_train,y_train,x_test,y_test):
tf.reset_default_graph()
sess = tf.InteractiveSession()
x = tf.placeholder(tf.float32, shape=[None, x_train.shape[1]])
y_ = tf.placeholder(tf.float32, shape=[None, 3])
W = tf.Variable(tf.zeros([x_train.shape[1],3]))
b = tf.Variable(tf.zeros([3]))
sess.run(tf.global_variables_initializer())
y = tf.matmul(x,W) + b
cross_entropy = tf.reduce_mean(
tf.nn.softmax_cross_entropy_with_logits(labels=y_, logits=y))
##train
train_step = tf.train.GradientDescentOptimizer(0.5).minimize(cross_entropy)
for i in range(1000):
train_step.run(feed_dict={x: x_train.tolist(), y_: onehot(y_train)})
#tf.argmax gives an index of the highest entry in a tensor along some axis
correct_prediction = tf.equal(tf.argmax(y,1), tf.argmax(y_,1))
#we can take this list of booleans and calculate the fraction correct
accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))
return accuracy.eval(feed_dict={x: x_test.tolist(), y_: onehot(y_test)})
sess.close()
In [60]:
def model_tf_nn(x_train,y_train,x_test,y_test):
tf.reset_default_graph()
n_inputs = x_train.shape[1]
n_hidden1 = 10
n_outputs = 3
learning_rate = 0.01
def neuron_layer(X, n_neurons, name, activation=None):
with tf.name_scope(name):
n_inputs = int(X.get_shape()[1])
stddev = 1 / np.sqrt(n_inputs)
init = tf.truncated_normal((n_inputs, n_neurons), stddev=stddev)
W = tf.Variable(init, name="weights")
b = tf.Variable(tf.zeros([n_neurons]), name="biases")
Z = tf.matmul(X, W) + b
if activation=="relu":
return tf.nn.relu(Z)
else:
return Z
X = tf.placeholder(tf.float32, shape=(None, n_inputs), name="X")
y = tf.placeholder(tf.int64, shape=(None), name="y")
with tf.name_scope("dnn"):
hidden1 = neuron_layer(X, n_hidden1, "hidden1", activation="relu")
logits = neuron_layer(hidden1, n_outputs, "output")
with tf.name_scope("loss"):
xentropy = tf.nn.sparse_softmax_cross_entropy_with_logits(labels=y, logits=logits)
loss = tf.reduce_mean(xentropy, name="loss")
with tf.name_scope("train"):
optimizer = tf.train.GradientDescentOptimizer(learning_rate)
training_op = optimizer.minimize(loss)
with tf.name_scope("eval"):
correct = tf.nn.in_top_k(logits, y, 1)
accuracy = tf.reduce_mean(tf.cast(correct, tf.float32))
init = tf.global_variables_initializer()
saver = tf.train.Saver()
n_epochs = 50
acc_test_high = 0.0
sess = tf.InteractiveSession()
sess.run(init)
for epoch in range(n_epochs):
for i in range(x_train.shape[0]):
x_data = x_train[i,:].reshape([1,x_train.shape[1]])
y_data = np.array([y_train[i]])
sess.run(training_op, feed_dict={X: x_data, y: y_data })
acc_train = accuracy.eval(feed_dict={X: x_train, y: y_train})
acc_test = accuracy.eval(feed_dict={X: x_test, y: y_test})
if acc_test > acc_test_high:
acc_test_high = acc_test
# print(epoch, "Train accuracy:", acc_train, "Test accuracy:", acc_test)
# save_path = saver.save(sess, "./one_layer.ckpt")
return acc_test_high
sess.close()
In [61]:
#data dir test
data_dir = '../data/greco_mri'
behav_dir = '../data/greco_behav'
#rois
rois = ['left_CA1','right_CA1','left_DG','right_DG'];
#subjec list
subjects = ['S1_A','S2_B','S3_A','S4_A','S5_A','S6_A','S7_A','S8_B',
'S9_A','S10_B','S11_B','S12_B','S13_B','S14_B','S15_B',
'S16_A','S21_B','S22_B','S24_A'];
subjects = np.array(subjects);
In [62]:
scores_svm = np.empty((subjects.shape[0],4))
scores_reg = np.empty((subjects.shape[0],4))
scores_nn = np.empty((subjects.shape[0],4))
for i,subj in enumerate(subjects):
for run in [0,1,2,3]:
behav = pd.read_table(join(behav_dir,subj,subj + '.txt'))
fname = join(data_dir,subj + '_right_DG.csv')
betas= pd.read_csv(fname,header=None)
test_ind = behav['run_num'] == run
train_ind = behav['run_num'] != run
y_test = behav['currCity'][test_ind].as_matrix()
y_test_r = np.empty(y_test.shape)
y_test_r[y_test == 1] = 0
y_test_r[y_test == 2] = 1
y_test_r[y_test == 3] = 2
y_train = behav['currCity'][train_ind].as_matrix()
y_train_r = np.empty(y_train.shape)
y_train_r[y_train == 1] = 0
y_train_r[y_train == 2] = 1
y_train_r[y_train == 3] = 2
x_test = center_normalize(betas[test_ind].as_matrix())
x_train = center_normalize(betas[train_ind].as_matrix())
scores_svm[i,run]=model_svm(x_train,y_train,x_test,y_test)
scores_reg[i,run] = model_tf_regression(x_train,y_train,x_test,y_test)
scores_nn[i,run] = model_tf_nn(x_train,y_train_r,x_test,y_test_r)
In [72]:
print('Logistic regression (mean accuracy): ' + str(np.mean(scores_reg)))
print('SVM (mean accuracy): ' + str(np.mean(scores_svm)))
print('Neural network (mean accuracy): ' + str(np.mean(scores_nn)))
In [76]:
from scipy import stats
print('Logistic regression (1samp ttest):')
print(stats.ttest_1samp(np.mean(scores_reg,1),.33))
print
print('SVM (1samp ttest):')
print(stats.ttest_1samp(np.mean(scores_svm,1),.33))
print
print('Neural network (1samp ttest)')
print(stats.ttest_1samp(np.mean(scores_nn,1),.33))
In [ ]: