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

Question

我（相对）是 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 次试验，那么？任何帮助将不胜感激。

Answer 1

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

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 个脉冲序列（事件序列）后，我们可以利用 R 中的 barcode 包：

#install.packages("barcode")

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

如果我在 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?

plot

r

raster

poisson

ggplot2