# Homework 2

## 1. Properties of the Poisson Random Processes

#### a)

Each neuron is spiking independently according to a Poisson process with rate $\lambda_1$ and $\lambda_2$ respectively. Let the random variables $T_1 \sim exp(\lambda_1)$ and $T_2 \sim exp(\lambda_2)$ correspond to the inter-spike intervals for each respective neuron. Let $T_r$ be the ISI for the raw spike train recorded on the electrode.

Then, for $t \geq 0$,
$$$
P(T*r > t) = P(T_1 > t \text{ and } T_2 > t) = P(T_1 > t)P(T_2 > t) = e^{-\lambda_1 t} e^{-\lambda_2 t} = e^{-(\lambda_1 + \lambda_2)t}
$$$
So, the cdf is $F*{T*r}(t) = 1 - e^{-(\lambda_1 + \lambda_2)t}$ for $t \geq 0$ and $F*{T_r}(t) = 0$ otherwise (since both $T_1$ and $T_2$ take only positive values). This is precisely the cdf of an exponential random variable with parameter $\lambda_1 + \lambda_2$.

#### b)

As shown in part a, the ISI distribution is exponential with parameter $\lambda_1 + \lambda_2$. Thus, the raw spike train is a Poisson process with rate $\lambda_1 + \lambda_2$.

#### c)

Let $N_1$ be the number of times the first neuron spikes within the 1 ms interval, and $N_2$ be the number of times the second neuron spikes within the 1 ms interval. Note that $N_1 \sim Poisson(\lambda_1 * 1)$ and $N_1 \sim Poisson(\lambda_2 * 1)$ (assuming that each rate is given in terms of spikes/ms).
Then,
$$$
P(N_1 = 1 \text{ and } N_2 =1) = P(N_1 = 1)P(N_2 = 1) = \left(e^{-\lambda_1}\frac{\lambda_1}{1!}\right)\left(e^{-\lambda_2}\frac{\lambda_2}{1!}\right) = \lambda_1 \lambda_2 e^{-(\lambda_1 + \lambda_2)}
$$$

## 2. Bayesian Analysis of the Poisson Process

#### a)

Suppose we are given observations $D = \lbrace x_1, x_2, \ldots, x_N\rbrace$ of a Poisson random variable characterized by $\text{Poi}(x \mid \lambda) = e^{-\lambda}\frac{\lambda^x}{x!}$ for $\lambda > 0$ and $x \in \lbrace 0, 1, 2, ... \rbrace$.

The log-likelihood function in this case is
$$$
\ell(\lambda) = \log{\left(\prod*{i=1}^{N} e^{-\lambda}\frac{\lambda^{x_i}}{x_i!}\right)} = \sum*{i=1}^{N} \log{\left(e^{-\lambda}\frac{\lambda^{x*i}}{x_i!}\right)} = \sum*{i=1}^{N} (-\lambda + x*i \log{\lambda} - \log{x_i!})
$$$
The log-likelihood function reaches its maximum when
$$$
0 = \frac{d\ell(\lambda)}{d\lambda} = \sum*{i=1}^{N} \left(-1 + \frac{x*i}{\lambda}\right) = -N + \sum*{i=1}^{N} \left(\frac{x*i}{\lambda}\right) \implies \lambda*{MLE} = \frac{1}{N} \sum_{i=1}^{N} x_i
$$$

Therefore, the ML-estimate in this case is just the sample mean of the given observations $D$.

#### b)

Assume the conjugate prior $p(\lambda) = \text{Gamma}(\lambda \mid a, b) \propto \lambda^{a-1} e^{-\lambda b}$. Then, the posterior distribution is
$$$
p(\lambda \mid D) \propto p(D \mid \lambda) p(\lambda) = \left(\prod*{i=1}^{N} e^{-\lambda}\frac{\lambda^{x_i}}{x_i!}\right) \left(\lambda^{a-1}e^{-\lambda b}\right )
$$$
$$$
= \frac{e^{-N\lambda}\lambda^{\sum*{i=1}^{N} x*i}}{\prod*{i=1}^{N} x*i!}\left(\lambda^{a-1}e^{-\lambda b}\right ) \propto e^{-\lambda(N+b)} \lambda^{a - 1 + \sum*{i=1}^{N} x_i}
$$$

Thus, the posterior takes the form of a $\text{Gamma}(a + \sum_{i=1}^{N} x_i, N+b)$

#### c)

The mean of the posterior is given by $\frac{a + \sum*{i=1}^{N} x_i}{N+b}$. So, as $a \to 0$ and $b \to 0$, the mean of the posterior tends to the sample mean $\frac{\sum*{i=1}^{N} x_i}{N}$, which is identical to the MLE found in part a.

## 3. Point Process Model Validation

#### a) Model Comparison

The Poisson model with given parameter has a higher log-likelihood given the data than the Gaussian model, so we conclude that it better fits the data.

#### b)

After training a Gaussian and Poisson model on the training spike count data, the log-likelihood of the Poisson model is higher than the log-likelihood for the Gaussian model on the validation set. So, we conclude that the Poisson model is a better fit for the data.

### c)

The above plot is constructed by computing the average spike rate in 0.05 second bins. Based on the plot, the rate of spikes increases gradually (and starts to taper off) over the 1 second interval. The rate is clearly not time-invariant, so a Poisson process with a time-varying rate would be preferable.

### d)

Based on the peri-stimulus time histogram (PSTH) above, the rate of spikes increases approximately linearly over the 1 second interval. The rate is clearly not time-invariant, so a Poisson process with a close to linear time-varying rate would be preferable.

These values are expected based on the plot above and the inhomogeneous Poisson process model. We expect the mean number of spikes (and the variance) during a 1-second interval to be $\int_0^1 \lambda(s) ds$. Based on the plot above, the area under the curve of $\lambda(t)$ is slightly larger than 4, which is consistent with the sample mean and variance of the count distribution found above. However, the count distribution mean and variance alone do not capture the inhomogeneous nature of the spike rate; for this, the PSTH-type approach is needed to show the changing rate over time.

One physiological process which could explain the observed time-varying rate is a gradual increase of a stimulus over the interval. If a certain stimulus was increased in intensity over the interval, we would expect the rate to also increase gradually, as seen here.

## 4. Describing Real Data

### a)

This neuron seems to have a "preferred" direction at around 140 degrees (and equivalently, due to symmetry, at around 320 degrees). The bimodal shape of this graph is due to the 180 degree rotational symmetry of the presented stimulus. An orientation of $\theta$ degrees is equivalent to an orientation of $\theta + 180$ degrees.

### b)

### c)

These Fano indices do not correspond with the expected for a Poisson process (expect F=1). The difference from a Fano Factor of 1 is likely due to the effect of refractory periods and non-linear spike rates in real-world data. These factors cause the variance to be larger than the mean of the observed counts per stimulus presentation.

### d)

These distributions roughly resemble Gamma distributions with $n>1$. The ISI distribution expected under a homogeneous Poisson process model is exponential. However, due to the time-varying nature of the spike rate and the fact that refractory times would limit the number of extremely short ISI lengths, the ISI distribution for real-world data (as in this problem) is better modeled as a Gamma distribution.

### e)

As shown in the above plots, after rescaling the ISIs (assuming a piecewise constant conditional intensity function), the ISI distributions appear to more closely resemble exponential distributions, but not perfectly so (expected due to fact that our piecewise constant assumption is likely not completely accurate)

Based on the Kolmogorov-Smirnov goodness of fit test, we reject the null hypothesis that the rescaled ISI distributions follow an exponential distribution. According to this test, the rescaled spike times are not modeled well by a homegeneous Poisson process. Furthermore, this indicates that our model of the spiking rate as piecewise constant was flawed; in reality, the rate likely varies with time within each interval.