読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

ワイブル分布に従う確率変数の和の分布の最尤推定

※初公開時は R のコードが間違っていて、うまくパラメータが求まっていなかった。(修正:7/28)

f:id:abrahamcow:20150728171354p:plain

# R code
shp1 <- 5
shp2 <- 5
scl1 <- 2
scl2 <- 8
n <-1000
Z <- rweibull(n,shp1,scl1) + rweibull(n,shp2,scl2)

LL <- function(par){
  m1 <- par[1]
  eta1 <- par[2]
  m2 <- par[3]
  eta2 <- par[4]

  f_X <- function(x) dweibull(x,m1,eta1)    
  f_Y <- function(y) dweibull(y,m2,eta2)
  integrand <- function(x,z) f_Y(z-x)*f_X(x)  #たたみこみ
  f_Z <- function(z) integrate(function(x,z) f_Y(z-x)*f_X(x),0,Inf,z)$value
  
  L <- sapply(Z,f_Z)
  
  sum(log(L))
}

res <- optim(c(5,2,5,8),LL,control=list(fnscale=-1))

hist(Z,breaks="FD",freq=FALSE)
f_X0 <- function(x) dweibull(x,shp1,scl1)    
f_Y0 <- function(y) dweibull(y,shp2,scl2)
integrand0 <- function(x,z) f_Y0(z-x)*f_X0(x)
f_Z0 <- function(z) integrate(integrand0,0,Inf,z)$value
f_Z0 <- Vectorize(f_Z0)
curve(f_Z0(x),add=TRUE,lty=2,lwd=1.5)

f_X <- function(x) dweibull(x,res$par[1],res$par[2])    
f_Y <- function(y) dweibull(y,res$par[3],res$par[4])
integrand <- function(x,z) f_Y(z-x)*f_X(x)
f_Z <- function(z) integrate(integrand,0,Inf,z)$value
f_Z <- Vectorize(f_Z)
curve(f_Z(x),add=TRUE,lwd=1.5)
legend("topleft",c("estimate","true"),lwd=1.5,lty=1:2)
> res
$par
[1] 5.911824 1.954939 5.021444 8.027303

$value
[1] -1954.873

$counts
function gradient 
     143       NA 

$convergence
[1] 0

$message
NULL

abrahamcow.hatenablog.com