2014-05-09 1 views
9

R에서 두 개의 확률 분포의 컨볼 루션을 계산하고 싶습니다. 도움이 필요합니다. 간단히하기 위해 mean = 1.0 및 stdev = 0.5로 정규 분포 된 변수 x와 평균 = 1.5 및 stdev = 0.75로 정규 분포 된 y를 가정 해 봅시다. 나는 z = x + y를 결정하고 싶다. 나는 z의 분포가 선험적으로 알려지지 않았 음을 이해합니다.convolution을 통해 두 개의 무작위 변수 추가하기 R

제가 작업하고있는 현실적인 예를 제외하고는 여러 배포본에 따라 배포되는 두 개의 임의의 변수에 추가가 필요합니다.

누구든지 x와 y의 확률 밀도 함수를 convuting하여 두 개의 임의 변수를 추가하는 방법을 알고 있습니까?

정상적으로 분포 된 임의의 값 (위의 매개 변수 사용)을 생성하고 로그에 정규 분포 된 임의의 값을 추가하려고했습니다. 그러나 대신 컨볼 루션 방법을 사용할 수 있는지 알고 싶습니다. 어떤 도움이라도 대단히 감사하겠습니다.

편집이 답변 주셔서 감사합니다. pdf를 정의하고 컨볼 루션 적분을 시도하지만 R은 통합 단계에서 불평합니다. f.t 기능이 경계하기 때문에 R 마지막 단계에서 불평

dlp3 <- function(x, a, b, g) { 
p1 <- 1/(x*abs(b) * gamma(a)) 
p2 <- ((log(x)-g)/b)^(a-1) 
p3 <- exp(-1* (log(x)-g)/b) 
d <- p1 * p2 * p3 
return(d) 
} 

f.m <- function(x) dlp3(x,3.2594,-0.18218,0.53441) 
f.s <- function(x) dlp3(x,9.5645,-0.07676,1.184) 

f.t <- function(z) integrate(function(x,z) f.s(z-x)*f.m(x),-Inf,Inf,z)$value 
f.t <- Vectorize(f.t) 
integrate(f.t, lower = 0, upper = 3.6) 

을 다음 내 통합 한도 아마 정확하지로 내 PDF 파일은 로그인 피어슨 3하고 있습니다. 이 문제를 해결하는 방법에 대한 아이디어가 있습니까?

+1

나는 당신이 [DISTR 패키지]를 체크 아웃 제안 (http://cran.r-project.org/web/packages/distr/index.html) 또는 적어도 빠른가에서 [비 네트 모양 받아 ] (http://cran.r-project.org/web/packages/distr/vignettes/newDistributions.pdf). 나는 그것이 당신이 찾고있는 것과 정확히 같다고 생각합니다. 무작위 값을 생성하는 데 사용한 전략이 완벽하게 유효하더라도. – MrFlick

답변

12

여기에는 한 가지 방법이 있습니다.

f.X <- function(x) dnorm(x,1,0.5)  # normal (mu=1.5, sigma=0.5) 
f.Y <- function(y) dlnorm(y,1.5, 0.75) # log-normal (mu=1.5, sigma=0.75) 
# convolution integral 
f.Z <- function(z) integrate(function(x,z) f.Y(z-x)*f.X(x),-Inf,Inf,z)$value 
f.Z <- Vectorize(f.Z)     # need to vectorize the resulting fn. 

set.seed(1)        # for reproducible example 
X <- rnorm(1000,1,0.5) 
Y <- rlnorm(1000,1.5,0.75) 
Z <- X + Y 
# compare the methods 
hist(Z,freq=F,breaks=50, xlim=c(0,30)) 
z <- seq(0,50,0.01) 
lines(z,f.Z(z),lty=2,col="red") 

패키지 distr를 사용

같은 것.

library(distr) 
N <- Norm(mean=1, sd=0.5)   # N is signature for normal dist 
L <- Lnorm(meanlog=1.5,sdlog=0.75) # same for log-normal 
conv <- convpow(L+N,1)    # object of class AbscontDistribution 
f.Z <- d(conv)     # distribution function 

hist(Z,freq=F,breaks=50, xlim=c(0,30)) 
z <- seq(0,50,0.01) 
lines(z,f.Z(z),lty=2,col="red") 
+0

dlp3 배포에 관해서는 방정식이 맞습니까? 당신의 함수'dlp3 (...)'이 갈라진다. 'x <- seq (0,100, .1);을 시도해보십시오. plot (x, dlp3 (x, 3, -0.2,0.5))'. 실제로, 큰 x에 대해 ~ x^4 * log (x)^2를 따라갈 수 있습니다. – jlhoward

관련 문제