hist.lockedTrain: Auto- and Cross-Intensity Function Estimate for Spike Trains
In STAR: Spike Train Analysis with R

Description Usage Arguments Details Value Note Author(s) References See Also Examples

hist.lockedTrain constructs and plot.hist.lockedTrain plots estimates of what Cox and Lewis (1966) call the auto- or cross-intensity functions. The auto-intensity function is also called the renewal density by Cox and Lewis and by Perkel et al (1967) while it is called the intensity of the Palm distribution by Ogata (1988). The (estimate of) the cross-intensity function is called cross-correlation function by Perkel et al (1967b) and cross-correlation histogram by Brillinger et al (1976).

## S3 method for class 'lockedTrain'
hist(x, bw, breaks = NULL, plot = TRUE, ...)
## S3 method for class 'hist.lockedTrain'
plot(x, style = c("Ogata", "Brillinger"),
                 CI, unit = "s", xlab, ylab, xlim, ylim,
                 type, pch, ...)

`x`	a `lockedTrain` object returned by the `lockedTrain` function.
`bw`	a non-negative numeric, the bin width.
`breaks`	a vector giving the breakpoints between cells. If `NULL` (default) breaks are built using argument `bw` and component `laglim` of `obj`.
`plot`	a logical. If `TRUE` a plot is generated as a side effect and nothing is returned, if `FALSE` a list of class `hist.lockedTime` is returned.
`style`	a character. The style of the plot, `"Ogata"` or `"Brillinger"`.
`CI`	a numeric vector with at most two elements. The coverage probability of the confidence intervals.
`unit`	a character. The unit in which the spike times are expressed.
`xlab`	a character. The x label. Default supplied.
`ylab`	a character. The y label. Default supplied.
`xlim`	a numeric. See `plot`. Default supplied.
`ylim`	a numeric. See `plot`. Default supplied.
`type`	see `lines`. Default supplied.
`pch`	see `plot`. Default supplied.
`...`	see `plot`.

The intensity of the Palm distribution (Ogata, 1988, p 13) is estimated by:

m(s) = Prob(event in (t+s,t+s+bw) | event at t) / bw

It is called renewal density by Perkel et al (1967) and defined by their Eq. 10, p 404. Under the null hypothesis of a stationary Poisson process it is a constant whose value is the mean discharge rate.

The cross-intensity function of two spike trains A and B is estimated by (Perkel et al, 1967b, p424, Eq. 4 and 5):

m(s;AB) = Prob(B event in (t+s,t+s+bw) | A event at t) / bw

The style argument of plot.hist.lockedTrain generates a plot looking like Fig. 6, p 18 of Ogata (1988) if set to "Ogata". Using style "Brillinger" plots the square root of the estimate.

When argument plot in hist.lockedTrain is set to FALSE a list of class hist.lockedTrain with the following components is returned:

`density`	the density estimate. Equivalent of the component `density` returned by `hist`.
`breaks`	a numeric vector with the breaks in between which spikes were counted. Similar to the component of the same name returned by `hist`.
`mids`	a numeric vector with the mid points of `breaks`. . Similar to the component of the same name returned by `hist`.
`bw`	the bin width used.
`nRef`	the total number of reference spikes used.
`refFreq`	the mean frequency of the reference neuron.
`testFreq`	the mean frequency of the test neuron.
`obsTime`	the total observation time used (in s).
`CCH`	a logical which is `TRUE` if a cross-intensity was estimated and `FALSE` in the case of an auto-intensity.
`call`	the matched call.

The confidence intervals are obtained from a Poisson distribution with parameter: refFreq * testFreq * bw * obsTime. Once the quantiles of the Poisson distribution have been obtained they are divided by: refFreq * bw * obsTime

These intervals are valid under the stationary Poisson null hypothesis for the auto-intensity estimates. They are valid under the stationary independent null hypothesis for the cross-intensity. There is NO NEED to assume that the test train is Poisson or renewal. See Perkel et al (1967b) and McFadden (1962) for a justification/proof of that. The square root transform of Brillinger (1976) and Brillinger et al (1976) is (in my opinion) a perfect example of shooting at a sparrow with a bazooka. An oversized method to get at an intuitively obvious result. There is moreover no need to stabilize the variance if your testing against a Poisson with a constant rate since then the variance of the null hypothesis is stable to start with. These (square root) transforms are useful for least square fits with a Poisson noise but NOT in the present context.

Christophe Pouzat christophe.pouzat@gmail.com

Ogata, Yosihiko (1988) Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes. Journal of the American Statistical Association 83: 9-27.

D. R. Cox and P. A. W. Lewis (1966) The Statistical Analysis of Series of Events. John Wiley and Sons.

J. A. McFadden (1962) On the Lengths of Intervals in a Stationary Point Process. Journal of the Royal Statistical Society. Series B, 24: 364-382

Perkel D. H., Gerstein, G. L. and Moore G. P. (1967) Neural Spike Trains and Stochastic Point Processes. I. The Single Spike Train. Biophys. J., 7: 391-418. http://www.pubmedcentral.nih.gov/articlerender.fcgi?tool=pubmed&pubmedid=4292791

Perkel D. H., Gerstein, G. L. and Moore G. P. (1967b) Neural Spike Trains and Stochastic Point Processes. I. Simultaneous Spike Trains. Biophys. J., 7: 419-440. http://www.pubmedcentral.nih.gov/articlerender.fcgi?tool=pubmed&pubmedid=4292792

David R. Brillinger, Hugh L. Bryant and Jose P. Segundo (1976) Identification of synaptic interactions. Biol Cybern, 22: 213-228.

David R. Brillinger (1976) Estimation of the Second-Order Intensities of a Bivariate Stationary Point Process. Journal of the Royal Statistical Society. Series B (Methodological), 38, 60-66.

varianceTime, renewalTestPlot, lockedTrain

## Reproduce Fig. 6 of Ogata 1988
data(ShallowShocks)
shalShocks <- lockedTrain(as.spikeTrain(ShallowShocks$Date),,c(0,500))
shalShocksH <- hist(shalShocks,5,plot=FALSE)
plot(shalShocksH,"Ogata",c(0.95,0.99),xlab="TIME LAG (DAYS)",ylab="NUMBER OF EVENTS PER DAY")
## Reproduce Fig. 7 of Ogata 1988
myBinNb <- 101
myMidPoints <- seq(from = 1, to = 6, length.out=myBinNb)
myMidPoints <- 10^myMidPoints/200
myBreaks <- c(0,myMidPoints[-myBinNb] + diff(myMidPoints) / 2)
shalShocksH2 <- hist(shalShocks,breaks=myBreaks,plot=FALSE)
yy <- abs(shalShocksH2$density - shalShocksH2$refFreq)
plot(shalShocksH2$mids[shalShocksH2$density>0],
     yy[shalShocksH2$density>0],
     pch = 1,
     xlim = c(0.001,10000),
     log = "xy",
     xlab = "TIME LAG (DAYS)",
     ylab = "NUMBER OF EVENTS PER DAY"
     )