#Se comienza estableciendo los valores de a, b y g0=P(S=0).
#Estos valores dependen de la distribución de la frecuencia.
#Se usará para tal fin la función abg0, donde
#dist = 'bin', 'binN' o 'poi' según la distribución de N.
#param = vector de parámetros de la distribución de N.
###########################################################
abg0 <- function(dist,param){
if(dist=='bin'){
p <- param[2]
n <- param[1]
a <- -p/(1-p)
b <- -(n+1)*a
g0 <- (1-p+p*f(0))^n
}
if(dist=='binN'){
p <- param[2]
n <- param[1]
a <- (1-p)
b <- (n-1)*a
g0 <- (p/(1-(1-p)*f(0)))^n
}
if(dist=='poi'){
a <- 0
b <- param
g0 <- exp(-param*(1-f(0)))
}
c(a,b,g0)
}
#La función panjer calcula la probabilidad siguiente:
# g_r = P(S=r), de forma recursiva.
# Debe estar definida previamente la función f(x)
# que corresponde a la función de probabilidad de las
# severidades.
###########################################################
panjer <- function(r,dist,param){
info <- abg0(dist,param)
a <- info[1]; b <- info[2]; g0 <- info[3]
g <- g0
C <- 1/(1-a*f(0))
for(i in 1:r){
aux <- 0
for(j in 1:i){
aux <- aux + C*(a+(i-j+1)*b/i)*f(i-j+1)*g[j]
}
g[i+1] <- aux
}
g[r+1]
}
#Ejemplo 1
###########################################################
#Severidades uniformes discretas con parámetro n = 10
#Frecuencia geométrica con parámetro p=1/4
f <- function(x){(1/10)*(x>=1 && x<=10)}
r <- 20
panjer(r,'binN',c(1,1/4))
## [1] 0.01557006
#Comprobación usando simulación
m <- 10000
i <- 0
uno_o_cero <- numeric(m)
set.seed(123)
N <- rgeom(m,1/4)
while(i<m){
i <- i + 1
if(N[i] == 0){
S <- 0
}else{
Y <- sample(1:10,N[i],replace = T)
S <- sum(Y)
}
uno_o_cero[i] <- 1*(S == r)
}
mean(uno_o_cero)
## [1] 0.0159
#Ejemplo 2
###########################################################
#Severidades Poisson con parámetro lambda = 7
#Frecuencia binomial con parámetros n= 10, p=1/3
f <- function(x){dpois(x,7)}
r <- 20
panjer(r,'bin',c(10,1/3))
## [1] 0.03408564
#Comprobación usando simulación
m <- 10000
i <- 0
uno_o_cero <- numeric(m)
set.seed(123)
N <- rbinom(m,10,1/3)
while(i<m){
i <- i + 1
if(N[i] == 0){
S <- 0
}else{
Y <- rpois(N[i],7)
S <- sum(Y)
}
uno_o_cero[i] <- 1*(S == r)
}
mean(uno_o_cero)
## [1] 0.034
#Ejemplo 3
###########################################################
#Severidades binomiales con parámetros n= 10, p=1/3
#Frecuencia Poisson con parámetro lambda = 7
f <- function(x){dbinom(x,10,1/3)}
r <- 20
panjer(r,'poi',7)
## [1] 0.04153745
#Comprobación usando simulación
m <- 10000
i <- 0
uno_o_cero <- numeric(m)
set.seed(123)
N <- rpois(m,7)
while(i<m){
i <- i + 1
if(N[i] == 0){
S <- 0
}else{
Y <- rbinom(N[i],10,1/3)
S <- sum(Y)
}
uno_o_cero[i] <- 1*(S == r)
}
mean(uno_o_cero)
## [1] 0.0404
No hay comentarios:
Publicar un comentario