This notebook presents how to perform maximum-likelihood parameter estimation for a single neuron.



In [1]:

    
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline



In [2]:

    
import sys
sys.path.append("../")
sys.path.append("../code/")



In [3]:

    
import filters
import likelihood_functions as lk
import PoissonProcessClasses as PP



In [4]:

    
import imp
imp.reload(filters)
imp.reload(lk)









    Out[4]:





<module 'likelihood_functions' from '../code/likelihood_functions.py'>



In [5]:

    
import sys
sys.path.append("../")
sys.path.append("../data/")

Reading input-output data:



In [6]:

    
# reading stimulus
Stim = pd.read_csv('../data/Stim.csv',header = None)
# reading location of spikes
tsp = np.hstack(np.array(pd.read_csv('../data/tsp.csv',header = None)))

Extracting a spike train from spike positions:



In [7]:

    
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)]).astype(int)

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}
H,H_orth,ht_domain = filters.createPostSpikeBasis(pars_h,dt)

Conditional Intensity (spike rate): $$\lambda_{\beta} = \exp(K(\beta_k)*Stim + H(\beta_h)*y)$$

Creating a matrix of covariates:



In [10]:

    
M = lk.construct_covariates(Stim,tsp,K,H_orth,dt)

Condtional intensity as a function of the covariates: $$ \lambda_{\beta} = \exp(M\beta) $$

Create a Poisson process model with this intensity:



In [11]:

    
model = PP.PPModel(M.T,dt = dt)

Setting initial parameters:



In [12]:

    
coeff_k0 = np.array([ 0.061453,0.284916,0.860335,1.256983,0.910615,0.488660,-0.887091,0.097441,0.026607,-0.090147])

coeff_h0 = np.zeros((5,))

pars0 = np.hstack((coeff_k0,coeff_h0))

Fitting the likelihood:



In [13]:

    
res = model.fit(y,start_coef = pars0)









    



../code/PoissonProcessClasses.py:167: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  if start_coef==None:
/Users/val/anaconda/lib/python3.4/site-packages/scipy/optimize/_minimize.py:362: RuntimeWarning: Method L-BFGS-B does not use Hessian information (hess).
  RuntimeWarning)



In [14]:

    
theta_MLE = res.x



In [ ]:



In [ ]: