import os import numpy as np import scipy as sp from scipy.linalg import sqrtm, inv import scipy.stats as stats import scipy.signal as signal import matplotlib.pyplot as plt import sympy import math %matplotlib inline def pupil_IRF(timepoints, s=1.0/(10**26), n=10.1, tmax=0.93): """ pupil_IRF defines the IRF (response to a transient input) of the pupil. Parameters ---------- timepoints: timepoints to evulate function s: scaling factor n: sets the width tmax: sets the time to peak Returns ------- y: IRF evaluated for 'timepoints' y_dt: IRF first derivative evaluated for 'timepoints' """ # in sympy: t = sympy.Symbol('t') y = ( (s) * (t**n) * (math.e**((-n*t)/tmax)) ) yprime = y.diff(t) # lambdify: y = sympy.lambdify(t, y, "numpy") y_dt = sympy.lambdify(t, yprime, "numpy") # evaluate: y = y(timepoints) y_dt = y_dt(timepoints) # normalize: y = y/np.std(y) y_dt = y_dt/np.std(y_dt) return (y, y_dt) # create the IRF: tmax_shifted = 0.93-0.3 sample_rate = 20 IRF_len = 3.0 # in seconds timepoints = np.linspace(0,IRF_len,IRF_len*sample_rate) IRF, IRF_prime = pupil_IRF(timepoints=timepoints) IRF_shifted, IRF_prime_shifted = pupil_IRF(timepoints=timepoints, tmax=tmax_shifted) # IRF_len = 20.0 # in seconds # timepoints = np.linspace(0,IRF_len,IRF_len*sample_rate) # IRF, IRF_prime = HRF(timepoints=timepoints) # plot the IRF: fig = plt.figure() plt.plot(timepoints, IRF_shifted, color='b') plt.plot(timepoints, IRF, color='r', ls='--') plt.plot(timepoints, IRF_prime, color='g', ls='--') plt.legend(['individual IRF', 'canonical IRF', 'canonical IRF dt',]) plt.title('Impulse Response Function') plt.xlabel('time (s)') plt.ylabel('a.u.') true_beta = 1 duration = 15 # in seconds times = np.array([2,2.5,5,6,8.75,9,10,10.5,]) input_signal = np.zeros(duration * sample_rate) for i in times: input_signal[i*sample_rate] = true_beta # convolve inputs with IRF: convolved_signal = (sp.convolve(input_signal, IRF, 'full'))[:-(IRF.shape[0]-1)] convolved_signal_shifted = (sp.convolve(input_signal, IRF_shifted, 'full'))[:-(IRF.shape[0]-1)] # # let's add some noise: # convolved_signal = convolved_signal + np.random.normal(0,0.2,len(convolved_signal)) # convolved_signal_shifted = convolved_signal_shifted + np.random.normal(0,0.2,len(convolved_signal_shifted)) # plot simulated convolved signal with noise: timepoints = np.linspace(0,duration,duration*sample_rate) fig = plt.figure() plt.plot(timepoints, convolved_signal_shifted, 'b') plt.plot(timepoints, convolved_signal, 'r', ls='--') plt.ylim(-1,8) plt.legend(['pupil time series', 'expected timeseries\nbased on canonical IRF',], loc=1) for i in times: plt.axvline(i, color='k', alpha=0.5) plt.title('pupil time series, with measurement noise') plt.xlabel('time (s)') plt.ylabel('a.u.') input_signal = np.zeros(duration * sample_rate) for i in times: input_signal[i*sample_rate] = 1 # convolve inputs with IRF: regr_convolved = (sp.convolve(input_signal, IRF, 'full'))[:-(IRF.shape[0]-1)] regr_convolved_dt = (sp.convolve(input_signal, IRF_prime, 'full'))[:-(IRF_prime.shape[0]-1)] # z-score measured data and regressors: regr_convolved = (regr_convolved - regr_convolved.mean()) / regr_convolved.std() regr_convolved_dt = (regr_convolved_dt - regr_convolved_dt.mean()) / regr_convolved_dt.std() convolved_signal_shifted = (convolved_signal_shifted - convolved_signal_shifted.mean()) / convolved_signal_shifted.std() fig = plt.figure() plt.plot(timepoints, regr_convolved, 'r',) plt.ylim(-4,6) plt.legend(['regressor'], loc=1) plt.title('model 1') plt.xlabel('time (s)') plt.ylabel('a.u.') fig = plt.figure() plt.plot(timepoints, regr_convolved, 'r',) plt.plot(timepoints, regr_convolved_dt, 'g') plt.ylim(-3,6) plt.legend(['regressor', 'regressor dt'], loc=1) plt.title('model 2') plt.xlabel('time (s)') plt.ylabel('a.u.') # make design matrix: designMatrix_0 = np.mat(regr_convolved).T designMatrix_1 = np.mat(np.vstack((regr_convolved, regr_convolved_dt))).T # multiple regression, yielding the betas: betas_0 = np.array(((designMatrix_0.T * designMatrix_0).I * designMatrix_0.T) * np.mat(convolved_signal_shifted).T).ravel() betas_1 = np.array(((designMatrix_1.T * designMatrix_1).I * designMatrix_1.T) * np.mat(convolved_signal_shifted).T).ravel() betas_1_combined = np.sign(betas_1[::2])*np.sqrt((betas_1[::2]**2) + (betas_1[1::2]**2)) # explained signal: explained_signal_0 = regr_convolved*betas_0[0] explained_signal_1 = regr_convolved*betas_1[0] + regr_convolved_dt*betas_1[1] # model fit: r_0, p_0 = sp.stats.pearsonr(convolved_signal_shifted, explained_signal_0) r_1, p_1 = sp.stats.pearsonr(convolved_signal_shifted, explained_signal_1) print 'model 0:' print '---------' print 'beta: {}\nR2: {}'.format(round(betas_0,3), round(r_0,3)) print print 'model 1:' print '---------' print 'beta: {} ({} ; {})\nR2: {}'.format(round(betas_1_combined,3), round(betas_1[0],3), round(betas_1[1],3), round(r_1,3)) print print '---------' print 'true beta: {}'.format(true_beta) print '---------' print # plotting: fig = plt.figure() plt.plot(timepoints, convolved_signal_shifted, 'b') plt.ylim(-3,6) plt.plot(timepoints, explained_signal_0, 'm', lw=1) plt.legend(['pupil timeseries', 'explained signal model 1'], loc=2) plt.title('predicted signal model 1') plt.xlabel('time (s)') plt.ylabel('a.u.') fig = plt.figure() plt.plot(timepoints, convolved_signal_shifted, 'b') plt.plot(timepoints, explained_signal_1, 'm', lw=1) plt.ylim(-3,6) plt.legend(['pupil timeseries', 'explained signal model 2'], loc=2) plt.title('predicted signal model 2') plt.xlabel('time (s)') plt.ylabel('a.u.') # residuals: residuals_0 = convolved_signal_shifted - explained_signal_0 residuals_1 = convolved_signal_shifted - explained_signal_1 fig = plt.figure() plt.plot(timepoints, residuals_0, 'c') plt.ylim(-3,6) plt.legend(['residuals model 1'], loc=1) plt.title('residuals') plt.xlabel('time (s)') plt.ylabel('a.u.') fig = plt.figure() plt.plot(timepoints, residuals_1, 'c') plt.ylim(-3,6) plt.legend(['residuals model 2'], loc=1) plt.title('residuals') plt.xlabel('time (s)') plt.ylabel('a.u.')