hist.lockedTrain: Auto- and Cross-Intensity Function Estimate for Spike Trains

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

View source: R/lockedTrain.R

Description

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).

Usage

1
2
3
4
5
6
## 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, ...)

Arguments

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.

Details

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.

Value

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.

Note

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.

Author(s)

Christophe Pouzat [email protected]

References

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.

See Also

varianceTime, renewalTestPlot, lockedTrain

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
## 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"
     )

STAR documentation built on May 31, 2017, 2:28 a.m.