2017-02-09 1 views
4

저는 r에서 화합물 포아송 과정을 시뮬레이션하려고합니다. 이 과정은 $ \ sum_ {j = 1}^{N_t} Y_j $에 의해 정의됩니다. $ Y_n $은 i.i.d와 독립적입니다. $ N (0,1) $ values와 $ N_t $는 $ 1 $ 매개 변수가있는 포아송 프로세스입니다. 난 행운을 빌어 r에서 이것을 시뮬레이트하려하고있다. Simutale 0부터 T까지 CPP :r에서 화합물 포아송 과정을 시뮬레이트합니다

초급자 : I는 다음과 같이 계산하는 알고리즘이 $ K = 0 $

반복하면서 $ \ sum_ {I = 1}^K T_i < T $으로

세트 $에서 K = K + 1 $

시뮬레이션 $ T_k \ 심 EXP (\ 람다) $ (내 케이스 $의 \ 람다 = 1 $)

시뮬레이션 $ Y_k \ SIM N (0 , 1) $ (이것은 특별한 경우 일 뿐이므로 이것을 배포본으로 변경할 수 있습니다.)

궤적은 다음과 같이 주어진다. 궤적은 다음과 같이 주어진다. 여기서, $ N (t) = sup (k : \ sum_ {i = 1}^k T_i \ leq t) $

프로세스를 계획 할 수 있도록 누군가를 r에서 시뮬레이션하도록 도와 줄 수 있습니까? 나는 시도했으나 끝내지 못했습니다.

+0

'rnorm','rexp '및'while'을 사용할 수 있습니까? 속도는 느리지 만 다른 프로그래밍 언어와는 전혀 다른 것입니다. 너 뭐 해봤 니? – AdamO

답변

2

N_t와 X_t를 결정하는 누적 합계에 cumsum을 사용하십시오. 이 예시 코드는 시뮬레이션 할 횟수를 지정하고 n, 시간을 n.tx의 값으로 시뮬레이트하고 (수행 한 내용을 표시하기 위해) 궤도를 그립니다.

n <- 1e2 
n.t <- cumsum(rexp(n)) 
x <- c(0,cumsum(rnorm(n))) 
plot(stepfun(n.t, x), xlab="t", ylab="X") 

Figure

이 알고리즘은 낮은 수준의 최적화 기능에 의존하기 때문에이 빠르고 : 나는 그것을 테스트 여섯 살짜리 시스템은 3 백만 이상 (시간, 값) 쌍을 생성합니다 초당. 시뮬레이션을위한 충분한 일반적으로 좋은,하지만 확실히 우리는 위의 코드를 활용할 수 시간 T. 밖으로 시뮬레이션을 생성하기 위해 묻는 문제, 만족하지 않지만,이 솔루션은 조금 까다 롭습니다


. 그것은 시간 T 이전에 포아송 프로세스에서 몇 번 발생할 것인지에 대한 합리적인 상한을 계산합니다. 도착 시간을 생성합니다. 이것은 시간 T에 실제로 도달하지 않은 (드문) 이벤트에서 프로 시저를 반복하는 루프에 랩핑됩니다.

추가 복잡성은 점근 계산 시간을 변경하지 않습니다.

T <- 1e2   # Specify the end time 
T.max <- 0   # Last time encountered 
n.t <- numeric(0) # Inter-arrival times 
while (T.max < T) { 
    # 
    # Estimate how many random values to generate before exceeding T. 
    # 
    T.remaining <- T - T.max 
    n <- ceiling(T.remaining + 3*sqrt(T.remaining)) 
    # 
    # Continue the Poisson process. 
    # 
    n.new <- rexp(n) 
    n.t <- c(n.t, n.new) 
    T.max <- T.max + sum(n.new) 
} 
# 
# Sum the inter-arrival times and cut them off after time T. 
# 
n.t <- cumsum(n.t) 
n.t <- n.t[n.t <= T] 
# 
# Generate the iid random values and accumulate their sums. 
# 
x <- c(0,cumsum(rnorm(length(n.t)))) 
# 
# Display the result. 
# 
plot(stepfun(n.t, x), xlab="t", ylab="X", sub=paste("n =", length(n.t))) 
+0

고마워요, 이건 엄청난 도움이됩니다. 더 좋은 줄거리를 얻을 수 있습니까? 즉 모든 점 주위에 원을 표시하지 않는 점은 무엇입니까? – Slangers

+0

예 :'? plot.stepfun'에있는 도움말 페이지를 참조하십시오. 그것은 마지막 줄에서'plot'을 호출 할 때'do.points = FALSE'라는 인자를 포함하도록 알려주고 플롯의 모양을 제어하기위한 다른 옵션을 설명합니다. – whuber

+0

대단히 감사합니다! :) – Slangers