我正在尝试为马尔可夫链生成一系列等待时间,其中等待时间是指数分布的数字,速率等于 1。但是,我不知道该过程的转换次数,而是该过程花费的总时间。
因此,例如:
t <- rexp(100,1)
tt <- cumsum(c(0,t))
t
是连续且独立的等待时间tt
的向量,是从 0 开始的实际转换时间的向量。
同样,问题是我不知道 t 的长度(即转换的数量),而不知道总等待时间将经过多少(即 tt 中最后一个条目的下限)。
在 R 中生成它的有效方法是什么?
我正在尝试为马尔可夫链生成一系列等待时间,其中等待时间是指数分布的数字,速率等于 1。但是,我不知道该过程的转换次数,而是该过程花费的总时间。
因此,例如:
t <- rexp(100,1)
tt <- cumsum(c(0,t))
t
是连续且独立的等待时间tt
的向量,是从 0 开始的实际转换时间的向量。
同样,问题是我不知道 t 的长度(即转换的数量),而不知道总等待时间将经过多少(即 tt 中最后一个条目的下限)。
在 R 中生成它的有效方法是什么?
The Wikipedia entry for Poisson process has everything you need. The number of arrivals in the interval has a Poisson distribution, and once you know how many arrivals there are, the arrival times are uniformly distributed within the interval. Say, for instance, your interval is of length 15.
N <- rpois(1, lambda = 15)
arrives <- sort(runif(N, max = 15))
waits <- c(arrives[1], diff(arrives))
Here, arrives
corresponds to your tt
and waits
corresponds to your t
(by the way, it's not a good idea to name a vector t
since t
is reserved for the transpose function). Of course, the last entry of waits
has been truncated, but you mentioned only knowing the floor of the last entry of tt
, anyway. If he's really needed you could replace him with an independent exponential (bigger than waits[N])
, if you like.
如果我做对了:您想知道填充您的时间间隔需要多少次转换。由于转换是随机且未知的,因此无法预测给定样本。以下是如何找到答案:
tfoo<-rexp(100,1)
max(which(cumsum(tfoo)<=10))
[1] 10
tfoo<-rexp(100,1) # do another trial
max(which(cumsum(tfoo)<=10))
[1] 14
现在,如果您预计需要绘制一些巨大的样本,例如rexp(1e10,1)
,那么也许您应该绘制“块”。抽取 1e9 个样本,看看是否sum(tfoo)
超过了您的时间阈值。如果是这样,请通过cumsum
. 如果没有,再抽取 1e9 个样本,以此类推。