In [1]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import random
%matplotlib inline
In [2]:
import sys
import os
sys.path.append(os.path.join(os.getcwd(),".."))
sys.path.append(os.path.join(os.getcwd(),"..","code"))
data_path = os.path.join(os.getcwd(),"..",'data')
sys.path.append(data_path)
In [3]:
import filters
import likelihood_functions as lk
import PoissonProcessClasses as PP
import auxiliary_functions as auxfun
In [4]:
# Reloading modules which are in development
import imp
imp.reload(filters)
imp.reload(auxfun)
imp.reload(lk)
imp.reload(PP)
Out[4]:
Reading input-output data:
In [5]:
# reading stimulus
Stim = np.array(pd.read_csv(os.path.join(data_path,'Stim.csv'),header = None))
# reading location of spikes
tsp = np.hstack(np.array(pd.read_csv(os.path.join(data_path,'tsp.csv'),header = None)))
Extracting a spike train from spike positions:
In [6]:
dt = 0.01
tsp_int = np.ceil((tsp - dt*0.001)/dt)
tsp_int = np.reshape(tsp_int,(tsp_int.shape[0],1))
tsp_int = tsp_int.astype(int)
y = np.array([item in tsp_int for item in np.arange(Stim.shape[0]/dt)+1]).astype(int)
Displaying a subset of the spike train:
In [7]:
fig, ax = plt.subplots(figsize=(16, 2))
fig = ax.matshow(np.reshape(y[:1000],(1,len(y[:1000]))),cmap = 'Greys',aspect = 15)
Creating filters:
In [8]:
# create a stimulus filter
kpeaks = np.array([0,round(20/3)])
pars_k = {'neye':5,'n':5,'kpeaks':kpeaks,'b':3}
K,K_orth,kt_domain = filters.createStimulusBasis(pars_k, nkt = 20)
In [9]:
# create a post-spike filter
hpeaks = np.array([0.1,2])
pars_h = {'n':5,'hpeaks':hpeaks,'b':.4,'absref':0.}
H,H_orth,ht_domain = filters.createPostSpikeBasis(pars_h,dt)
In [10]:
# Interpolate Post Spike Filter
MSP = auxfun.makeInterpMatrix(len(ht_domain),1)
MSP[0,0] = 0
H_orth = np.dot(MSP,H_orth)
In [11]:
MSP.shape
Out[11]:
Conditional Intensity (spike rate):
$$\lambda_{\beta} = \exp(K(\beta_k)*Stim + H(\beta_h)*y + dc)$$($\beta_k$ and $\beta_h$ are the unknown coefficients of the filters and $dc$ is the direct current).
Since the convolution is a linear operation the intensity can be written in the following form:
$$\lambda_{\beta} = \exp(M_k \beta_k + M_h\beta_h + \textbf{1}dc),$$where $M_k$ and $M_h$ are matrices depending on the stimulus and the response correspondingly and $\textbf{1}$ is a vector of ones.
Creating a matrix of covariates:
In [12]:
M_k = lk.construct_M_k(Stim,K,dt)
In [13]:
M_h = lk.construct_M_h(tsp,H_orth,dt,Stim)
Combining $M_k$, $M_h$ and $\textbf{1}$ into one covariate matrix:
In [14]:
M = np.hstack((M_k,M_h,np.ones((M_h.shape[0],1))))
The conditional intensity becomes: $$ \lambda_{\beta} = \exp(M\beta) $$
($\beta$ contains all the unknown parameters).
Create a Poisson process model with this intensity:
In [15]:
model = PP.PPModel(M.T,dt = dt/100)
Setting initial parameters for optimization:
In [16]:
coeff_k0 = np.array([-0.02304,
0.12903,
0.35945,
0.39631,
0.27189,
0.22003,
-0.17457,
0.00482,
-0.09811,
0.04823])
coeff_h0 = np.zeros((5,))
dc0 = 0
pars0 = np.hstack((coeff_k0,coeff_h0,dc0))
In [17]:
# pars0 = np.hstack((np.zeros((10,)),np.ones((5,)),0))
Fitting the likelihood (here using Limited Memory BFGS method with 500 iterations):
In [18]:
res = model.fit(y,start_coef = pars0,maxiter = 500,method = 'L-BFGS-B')
Optimization results:
In [19]:
print(res)
Creating the predicted filters:
In [20]:
k_coeff_predicted = res.x[:10]
h_coeff_predicted = res.x[10:15]
In [21]:
kfilter_predicted = np.dot(K,k_coeff_predicted)
hfilter_predicted = np.dot(H_orth,h_coeff_predicted)
In [22]:
k_coeff = np.array([ 0.061453,0.284916,0.860335,1.256983,0.910615,0.488660,-0.887091,0.097441,0.026607,-0.090147])
h_coeff = np.array([-15.18,38.24,-67.58,-14.06,-3.36])
kfilter_true = np.dot(K,k_coeff)
hfilter_true = np.dot(H_orth,h_coeff)
In [23]:
plt.plot(-kt_domain[::-1],kfilter_predicted,color = "r",label = 'predicted')
plt.hold(True)
plt.plot(-kt_domain[::-1],kfilter_true,color= "blue",label = 'true')
plt.hold(True)
plt.plot(-kt_domain[::-1],np.dot(K,coeff_k0),color = "g",label = 'initial')
plt.legend(loc = 'upper left')
Out[23]:
In [ ]:
In [24]:
plt.plot(ht_domain,hfilter_predicted,color = "r",label = 'predicted')
plt.hold(True)
plt.plot(ht_domain,hfilter_true,color = "b",label = 'true')
plt.hold(True)
plt.plot(ht_domain,np.dot(H_orth,coeff_h0),color = "g",label = 'initial')
plt.legend(loc = 'lower right')
Out[24]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: