如果我在 R 中的总持续时间内有一系列尖峰时间，我如何获得非均匀尖峰的光栅图？答案

【问题标题】：How do I get a raster plot of nonhomogeneous spikes if I have a sequence of spike times over the total duration, in R?如果我在 R 中的总持续时间内有一系列尖峰时间，我如何获得非均匀尖峰的光栅图？
【发布时间】：2018-04-25 12:40:18
【问题描述】：

我是 R 新手（相对而言），所以这个问题可能有一个简单的答案。我一直在尝试在 R 中绘制非齐次泊松过程。通过下面的代码，在一次 10 秒的试验中，我有一个序列nhpp1，它由在这个特定试验中出现尖峰的时间戳组成。我如何获取这个序列nhpp1 并获得它的光栅图。更重要的是，我想重复（复制？）整个事情进行 10 次试验，并得到一个看起来像这样的栅格图：（请看下面的代码）

 nhpp <- function(lambda){
 set.seed(1)
 t_max = 10
 t = 0 
 s = 0
 Lambda <- function(tupper) integrate(f=lambda, lower =0, upper= 
           tupper)$value
 LambdaInv <- function(s) {
         v <- seq(0, t_max+1, length=1000)
         min(v[which(Vectorize(Lambda)(v) >= s)])
           }
 X = numeric(0)
 while(t <= t_max){
  u <- runif(1)
  s <- s-log(u)
  t <- LambdaInv(s)
  X <- c(X,t)

  }
 return(X)
}

lambda <- function(t) 100*(sin(pi*t)+1)
nhpp1 <- nhpp(lambda)

我已经有了尖峰时间戳，我需要帮助来找到一种方法来绘制这个试验（在时间轴上出现尖峰的地方出现小条）以及如何复制这个过程进行 10 次试验？任何帮助将不胜感激。

【问题讨论】：

标签： r plot ggplot2 raster poisson

【解决方案1】：

我使用另一个基于 Poisson 模型的脉冲序列生成器。

arrival_time_v3 <- function(firing_rate,tMax,sampling_rate){
lambda <- firing_rate/sampling_rate

## find the number 'n' of exponential r.vs required by imposing that
## Pr{N(t) <= n} <= 1 - eps for a small 'eps'
n <- qpois(1 - 1e-8, lambda = lambda * tMax)

## simulate exponential interarrivals the
X <- rexp(n = 2*n, rate = lambda)
S <- c(0, cumsum(X))

arr_time <- S[which(S <= tMax)]
arr_time <- as.integer(arr_time)
arr_time <- arr_time[which(arr_time!=0)]
arr_time <- unique(arr_time)
return(arr_time)
}

num_of_spike_trains <- 10
firing_rate_arr <- matrix(c(10,20,30,40,50,60,70,80,90,100),1,num_of_spike_trains)
sampling_rate <- 10000
durartion_sample <- sampling_rate*10 ##10 sec

spike_arrival_time_mat_list <- list()
for(i in seq(1,num_of_spike_trains,1)){
spike_arrival_time_mat_list[[i]] <- t(as.matrix(arrival_time_v3(firing_rate_arr[i],durartion_sample,sampling_rate)))
}

生成10个spike trains（events trains）后，我们可以在R中使用barcode包：

#install.packages("barcode")

library(barcode)
barcode(spike_ariv_time_mat_list,xlab="Time",main="Spike Trains")

【讨论】：