Как получить растровую диаграмму неоднородных всплесков, если у меня есть последовательность времен всплесков по всей продолжительности в R? - PullRequest
Новая
       32

Как получить растровую диаграмму неоднородных всплесков, если у меня есть последовательность времен всплесков по всей продолжительности в R?

0 голосов
/ 25 апреля 2018

Я новичок в 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)

Raster Plot

У меня уже есть метки времени всплеска, мне нужна помощь в поиске способа построения этого одного испытания (с крошечными столбиками, возникающими там, где всплески происходят на временной шкале) и как повторить этот процесс для 10 испытаний, затем?Любая помощь будет очень ценится.

1 Ответ

0 голосов
/ 27 сентября 2018

Я использую другой генератор поездов с шипами, основанный на модели Пуассона. 

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 поездов с шипами (событийных поездов) мы можем использовать пакет barcode в R: 

#install.packages("barcode")

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

enter image description here

...