Modelo Beta Prime

Introdución

La distribución beta prime (BP) es una distribución de dos parámetros en la línea real positiva, que puede ser interpretada como la distribución de la razón de posibilidades de una variable distribuida según la distribución beta, es decir, si \(X \sim Beta(\alpha, \beta)\), entonces \(Y = \displaystyle \frac{X}{1-X} \sim BP(\alpha, \beta)\) con \(\alpha > 0\) y \(\beta >0\) estamos adoptando la parametrización de la distribución BP en términos de la media y los parámetros de precisión que fue propuesta por Bourguignon et al. (2018).

Una ventaja de usar esta parametrización es que podemos introducir estructuras de regresión para cada media y parámetros de precisión y la interpretación de los coeficientes de regresión es sencilla en términos de ellos como en los modelos lineales generalizados.

Distribución

La variable aleatoria \(Y\) se distribuye BP (Bourguignon et al. (2018)) con función de densidad de probabilidad dada por: \[ f(y;\mu , \phi)= \frac{y^{\mu(\phi+1)-1}(1+y)^{[-\mu(\phi+1)+\phi+2]}}{B(\mu(1+\phi),\phi+2)},y>0\]

donde \(\mu > 0\) y \(\phi > 0\) son parámetros de media y precisión, respectivamente, \(B(\mu(1+\phi),\phi+2) = \Gamma (\mu(1+\phi))\Gamma(\phi+2)\) es la función beta y \(\gamma(\mu(1+\phi)) = \int_{0}^{\infty} \omega^{\mu(1+\phi)-1}e^{-\omega} \cdot d\omega\) es la función gamma. De ahora en adelante, usamos la \(Y \sim BP(\mu,\phi)\) para indicar que \(Y\) es un variable aleatoria que sigue una distribución BP. La media y la varianza de \(Y\) son: \[ E[Y]=\mu , Var[Y]=\frac{\mu(1+\mu)}{\phi}\]

Modelo de regresión Beta prime

Considere \(n\) variables aleatorias independientes \(Y_{1}, . . . , Y_n\) donde cada \(Y_i\) , \(i = 1, . . . , n\) tiene PB distribución con media \(\mu_i\) y parámetro de precisión \(\phi_i\).\ Bourguignon et al. (2018) propusieron el modelo de regresión BP que se define por dos relaciones funcionales \[g_1(\mu_i)= \eta_{1i}=x_i^T \beta ,g_2(\mu_i)= \eta_{2i}=z_i^T v\]

Donde \(g_1:\mathbb{R}\longrightarrow \mathbb{R}^+\) y \(g_2:\mathbb{R}\longrightarrow \mathbb{R}^+\) son funciones de enlace estrictamente monótonas, positivas y al menos dos veces diferenciables, \(\eta_{1i}\) y \(\eta_{2i}\) son los predictores lineales, \(\beta=(\beta_1,...\beta_p)^T(\beta \in \mathbb{R}^p)\) y \(V = (v_1, ..., v_q)^T (V \in \mathbb{R}^q,q<n-p)\) son vectores de parámetros desconocidos a estimar, y $x_i=(x_{i1}, … ,x_{ip})^T $ y \(z_i=(z_{i1}, ... ,z_{iq})^T\) son observaciones sobre p y q regresores conocidos, para \(i = 1, . . . , n\). Además, suponemos que las matrices de covariables \(X = (x_1 ,...,x_n)^T\) y \(z = (z_1 ,...,z_n)^T\) tienen rango p y q, respectivamente.

MLE

En la literatura existen pocos trabajos que traten de la distribución de la PB. McDonald (1987) discutió sus propiedades y obtuvo las estimaciones de máxima verosimilitud (ML) de los parámetros del modelo. Las versiones corregidas de sesgo de los MLE de los parámetros que indexan la distribución de BP fueron obtenido por \(Stosíc\) y \(Cordeiro (2009)\). Vale la pena mencionar que todos los trabajos relacionados anteriormente tienennconsiderado la parametrización habitual de la distribución de la BP. Teniendo en cuenta la parametrización adoptamos, Bourguignon et al. (2018) utilizaron el método ML para estimar los parámetros que indexan el modelo de regresión de BP. Sin embargo, en muestras de pequeño tamaño, los estimadores ML de estos parámetros (especialmente para estructura de precisión) puede ser extremadamente parcial. Por lo tanto, es importante considerar estimadores alternativos con sesgos más pequeños.

La función logarítmica de verosimilitud para \((\beta, V)\) dados los valores observados \(y_1, . . . , y_n\) es \[\ell(\beta,V)=\sum_{i=1}^{n}\ell(\mu_i, \phi_i) \] siendo \[ \ell(\mu_i,\phi_i)= [\mu_i(1+\phi_i)-1]log(y_i)-[\mu_i(1+\phi_i)+\phi_i + 2]log(1+y_i)\] \[log[\Gamma(\mu_i(1+\phi_i))]-log[\Gamma(\phi_i+2)]+log[\Gamma(\mu_i(1+\phi_i)+\phi_i+2)]\] \(\mu_i=g_1^{-1}(\eta_{1i})\) y \(\mu_i=g_2^{-1}(\eta_{2i})\) son funciones de \(\beta\) y \(V\), respectivamente como se define en el modelo de regresión. Se describe un método para obtener las estimaciones de los parámetros del modelo BP definido por la distribución y en la regresión del modelo se describe en detalle en Bourguignon et al. (2018)

MLE en R

####################################################################
#
# Librerias
#
###################################################################@

library(Formula)
library(LaplacesDemon)
library(ggplot2)
library(extraDistr)
library(gamlss)

library(MASS)
library(Matrix)
library(nleqslv)
############################################################################
#              Función de densidad de probabilidad
###########################################################################
dBP <- function(x,mu=1,sigma=1,log=FALSE)
{
  if (any(mu < 0)) 
    stop(paste("mu debe ser positivo", "\n", ""))
  if (any(sigma < 0)) 
    stop(paste("sigma debe ser positivo", "\n", ""))
  if (any(x <= 0)) 
    stop(paste("x debe ser positivo", "\n", ""))
  
  a <- mu*(1+sigma)
  b <- 2 + sigma 
  
  fy <- dbetapr(x, shape1 = a, shape2 = b, scale = 1, log = log)
  fy
  
}
#############################################################################
#               Función de distribución acumulativa 
############################################################################
pBP <-  function(q,mu=1,sigma=1, lower.tail = TRUE, log.p = FALSE)
{
  if (any(mu < 0)) 
    stop(paste("mu debe ser positivoe", "\n", ""))
  if (any(sigma < 0)) 
    stop(paste("sigma debe ser positivo", "\n", ""))
  if (any(q < 0)) 
    stop(paste("q debe ser positivo", "\n", ""))
  
  a <- mu*(1+sigma)
  b <- 2 + sigma
  
  cdf <- pbetapr(q, shape1 = a, shape2 = b, scale=1, lower.tail = lower.tail, 
               log.p = log.p)
  cdf
}
#############################################################################
#                              Números Aleatorios
#############################################################################
rBP <- function(n,mu=1,sigma=1)
{
  if (any(mu < 0)) 
    stop(paste("mu debe ser positivo", "\n", ""))
  if (any(sigma < 0)) 
    stop(paste("sigma debe ser positivo", "\n", ""))
  if (any(n <= 0)) 
    stop(paste("n debe ser entero positivo", "\n", ""))
  
  n <- ceiling(n)
  
  a <- mu*(1+sigma)
  b <- 2 + sigma  
  
  r <- rbetapr(n,shape1=a,shape2=b,scale=1)
  
  r
}
#############################################################################
#                           Función Cuantil
############################################################################
qBP <- function(p, mu = 1, sigma = 1, lower.tail = TRUE, log.p = FALSE)
{
  if (any(mu < 0)) 
    stop(paste("mu debe ser positivo", "\n", ""))
  if (any(sigma < 0)) 
    stop(paste("sigma debe ser positivo", "\n", ""))
  if (any(n <= 0)) 
    if (any(p <= 0) | any(p >= 1)) 
      stop(paste("p debe estar entre 0 y 1", "\n", ""))  
  
  a <- mu*(1+sigma)
  b <- 2 + sigma 
  
  q <- qbetapr(p, shape1 = a, shape2 = b,scale=1, lower.tail = lower.tail, log.p = log.p)
  
  q
} 
#############################################################################
#                           familia
############################################################################
BP <- function (mu.link = "log", sigma.link = "log") 
{
  mstats <- checklink("mu.link", "Beta Prime", substitute(mu.link), c("log", "identity", "sqrt"))
  dstats <- checklink("sigma.link", "Beta Prime",substitute(sigma.link), c("log", "identity", "sqrt"))
  structure(list(family = c("BP", "Beta Prime"), 
                 parameters = list(mu = TRUE,sigma = TRUE), nopar = 2, type = "Continuous", 
                 mu.link = as.character(substitute(mu.link)), 
                 sigma.link = as.character(substitute(sigma.link)), 
                 mu.linkfun = mstats$linkfun, 
                 sigma.linkfun = dstats$linkfun, 
                 mu.linkinv = mstats$linkinv, 
                 sigma.linkinv = dstats$linkinv, 
                 mu.dr = mstats$mu.eta, 
                 sigma.dr = dstats$mu.eta, 
                 
                 
                 dldm = function(y, mu, sigma){
                 a <- mu*(1+sigma)
                 b <- mu*(1+sigma)+sigma+2
                 Phi <-  (1+sigma)  
                 yast <- log(y) - log(1+y) 
                 muast <- digamma(a) - digamma(b) 
                 dldm <- Phi*(yast - muast)
                 
                 dldm
                 }, 
                 d2ldm2 = function(mu, sigma){ 
                   Phi2 <- (1+sigma)^2
                   a <- mu*(1+sigma)
                   b <- mu*(1+sigma)+sigma+2
                   d2dldm2 <- -Phi2*(trigamma(a) - trigamma(b))
                   
                   d2dldm2
                   },
                 dldd = function(y, mu, sigma){ 
                   Phi <-  (1+sigma)
                   a <- mu*(1+sigma)
                   b <- mu*(1+sigma)+sigma+2 
                   ystar <- mu*log(y) - (1+mu)*log(1+y) 
                   mustar <- mu*digamma(a) - (1+mu)*digamma(b) + digamma(Phi+1)
                   
                   dldd <- ystar - mustar
                   
                   dldd
                  }, 
                 d2ldd2 = function(mu,sigma){
                   Phi <-  (1+sigma)
                   a <- mu*(1+sigma)
                   b <- mu*(1+sigma)+sigma+2 
                   
                   d2ldd2 <- -(mu^2)*trigamma(a) + ((1+mu)^2)*trigamma(b) - trigamma(Phi+1)
                   
                   d2ldd2
                   
                }, 
                 d2ldmdd = function(mu,sigma){
                   
                   a <- mu*(1+sigma)
                   b <- mu*(1+sigma)+sigma+2
                   Phi <-  (1+sigma)  
                   gammaast <- Phi*(trigamma(b) + mu*(trigamma(b)-trigamma(a)))
                   
                   d2ldmdd <- gammaast
                   
                   d2ldmdd
                
                   }, 
                 G.dev.incr = function(y, mu, sigma,...){-2*dBP(y, mu, sigma, log = TRUE)}, 
                 rqres = expression(rqres(pfun = "pBP", type = "Continuous", y = y, mu = mu, sigma = sigma)), 
                 mu.initial = expression({mu <- mean(y)}), 
                 sigma.initial = expression({sigma <-  mean(y)*(1+mean(y))/var(y) }), 
                 mu.valid = function(mu) all(mu > 0), 
                 sigma.valid = function(sigma) all(sigma > 0), 
                 y.valid = function(y) all(y > 0)),
                class = c("gamlss.family","family"))
}
#####################################################################
#          BIAS OF THE BETA AND NU - COX-SNELL
#####################################################################
bias <- function(x){ 
  beta1 = matrix(x[1:p], ncol=1)
  nu1 = matrix(x[-(1:p)], ncol=1)
  eta1.est = X%*%beta1
  eta2.est = Z%*%nu1
      
    mu.est <- as.vector(exp(eta1.est))
    phi.est <- as.vector(exp(eta2.est))

  dg1 <- psigamma(mu.est*(1+phi.est), deriv=1)
    dg11 <- psigamma(mu.est*(1+phi.est) + phi.est + 2, deriv=1)
  dg12 <- psigamma(phi.est + 2, deriv=1)
    dg2 <- psigamma(mu.est*(1+phi.est), deriv=2)
    dg21 <- psigamma(mu.est*(1+phi.est) + phi.est + 2, deriv=2)
  dg22 <- psigamma(phi.est + 2, deriv=2)

    a1 <- dg1 - dg11
  b1 <- mu.est^2*dg1 - (1+mu.est)^2*dg11 + dg12
    c1 <- dg2 - dg21
    d1 <- (1+mu.est)^2*dg21 - mu.est^2*dg2
    e1 <- (1+mu.est)^3*dg21 - mu.est^3*dg2 - dg22

  m1 <- -( (1+phi.est)^3*c1*(mu.est)^3 )/2 - ( (1+phi.est)^2*a1*mu.est*mu.est )/2 
    m2 <- ( (1+phi.est)*(-a1*mu.est + dg11)*mu.est*phi.est )/2 - ( (1+phi.est)^2*(mu.est*c1 - dg21)*mu.est^2*phi.est )/2
    m3 <- - ( (1+phi.est)*( 2*a1 + (1+phi.est)*mu.est*c1 - (1+phi.est)*dg21 )*mu.est^2*phi.est )/2  - ( (1+phi.est)*(dg11 - a1*mu.est)*mu.est*phi.est )/2
    m4 <- ( (1+phi.est)*d1 + 2*dg11 - 2*mu.est*a1 )*mu.est*phi.est^2/2 - ( (1+phi.est)*(dg11 - mu.est*a1)*mu.est*phi.est )/2
    m5 <- ( (1+phi.est)*d1/2 )*phi.est^2*mu.est + ( (1+phi.est)*(dg11 - mu.est*a1)*mu.est*phi.est )/2
    m6 <- (e1*phi.est^3 - b1*phi.est*phi.est)/2

  M1 <- diag(m1)
    M2 <- diag(m2)
    M3 <- diag(m3)
    M4 <- diag(m4)
    M5 <- diag(m5)
    M6 <- diag(m6)
    
  E1 <- diag((1+phi.est)^2*(a1)*mu.est^2)
  E2 <- diag((1+phi.est)*(mu.est*a1 - dg11)*mu.est*phi.est)
  E4 <- diag(b1*phi.est^2)

  K11 <- t(X)%*%E1%*%X
    K12 <- t(X)%*%E2%*%Z
    K21 <- t(K12)
    K22 <- t(Z)%*%E4%*%Z
  k1 <- cbind(K11, K12)
  k2 <- cbind(K21, K22)
  IE <- rbind(k1, k2)
  K <- solve(IE)

    Kbb <- K[1:p,1:p]
  Kbnu <- K[1:p,(p+1):(p+q)]
  Knub <- t(Kbnu)
  Knunu <- K[(p+1):(p+q),(p+1):(p+q)]
    Pbb <- matrix(diag(X%*%Kbb%*%t(X)))
    Pbnu <- matrix(diag(X%*%Kbnu%*%t(Z)))
    #Pnub <- t(Pbnu)
    Pnunu <- matrix(diag(Z%*%Knunu%*%t(Z)))
    
  biasbeta <- Kbb%*%t(X)%*%(M1%*%Pbb +(M2+M3)%*%Pbnu + M5%*%Pnunu) + Kbnu%*%t(Z)%*%(M2%*%Pbb +(M4+M5)%*%Pbnu + M6%*%Pnunu)

    biasnu <- Knub%*%t(X)%*%(M1%*%Pbb +(M2+M3)%*%Pbnu + M5%*%Pnunu) + Knunu%*%t(Z)%*%(M2%*%Pbb +(M4+M5)%*%Pbnu + M6%*%Pnunu)
      
  return(c(biasbeta, biasnu))
}
#####################################################################
#          VECTOR SCORE - DAVID FIRTH
#####################################################################
score_DF <- function(x, v){
  U           = numeric(p+q)
  beta.df     = matrix(x[1:p], ncol=1)
  nu.df       = matrix(x[-(1:p)], ncol=1)
  eta1.est.df = X%*%beta.df
  eta2.est.df = Z%*%nu.df
      
    mu.est.df   = as.vector(exp(eta1.est.df))
    phi.est.df  = as.vector(exp(eta2.est.df))
   
  dg1         = psigamma(mu.est.df*(1+phi.est.df), deriv=1)
    dg11        = psigamma(mu.est.df*(1+phi.est.df) + phi.est.df + 2, deriv=1)
  dg12        = psigamma(phi.est.df + 2, deriv=1)

  a1 = dg1 - dg11
  b1 = mu.est.df^2*dg1 - (1+mu.est.df)^2*dg11 + dg12

    E1 = diag((1+phi.est.df)^2*(a1)*mu.est.df^2)
  E2 = diag((1+phi.est.df)*(mu.est.df*a1 - dg11)*mu.est.df*phi.est.df)
  E4 = diag(b1*phi.est.df^2)

  K11 = t(X)%*%E1%*%X
    K12 = t(X)%*%E2%*%Z
    K21 = t(K12)
    K22 = t(Z)%*%E4%*%Z
  k1  = cbind(K11, K12)
  k2  = cbind(K21, K22)
  IE  = rbind(k1, k2)
   
  A        = matrix(x, ncol=p+q)
  res.bias = bias(A)
  vies     = matrix(c(res.bias), ncol=1) 
    
  M       = IE%*%vies
  PHI     = diag(1+phi.est.df)
  D1      = diag(mu.est.df)
  y.star  = matrix(log(v/(1+v)), ncol=1)
  mu.star = matrix(psigamma(mu.est.df*(1+phi.est.df), deriv=0) - psigamma(mu.est.df*(1+phi.est.df) + phi.est.df + 2, deriv=0), ncol=1)
  D2      = diag(phi.est.df)
  v.nu    = matrix(mu.est.df*(y.star-mu.star) + psigamma(mu.est.df*(1+phi.est.df) + phi.est.df + 2, deriv=0) - log(1+v) - psigamma(2+phi.est.df, deriv=0), ncol=1 )
 
  U.beta  = t(X)%*%PHI%*%D1%*%(y.star-mu.star)
  U.nu    = t(Z)%*%D2%*%v.nu
  U       = c(as.vector(U.beta) - M[1:p] , as.vector(U.nu) - M[(p+1):(p+q)] )
    # print(U)
  U
    }
library(tidyverse)
library(lubridate)
library(ggthemes)
library(ggtext)
###############################################################################
#             Datos 
###############################################################################

library(tidyverse)
csv_data_BP <- read_csv('midata.xls')

Kilometres <-  csv_data_BP$Kilometres
Zone <- csv_data_BP$Zone
Bonus <- csv_data_BP$Bonus
Make <- csv_data_BP$Make
Insured <- csv_data_BP$Insured
Claims <- csv_data_BP$Claims
Payment <- csv_data_BP$Payment

rb = 1000              # Repeticiones para Bootstrap
p = 6                  # dimensiones de beta
q = 3                  # dimensiones de gamma
n  = length(Payment/Claims)
X  = matrix(c(rep(1,n), Kilometres, Zone, Bonus, Make, Insured), ncol = p)
Z  = matrix(c(rep(1,n), Kilometres, Zone), ncol = q)
Dobs = Payment/Claims
########################################################################################
#                         M.A.I.N   F.U.N.C.T.I.O.N
########################################################################################
conh = gamlss.control(trace = FALSE, autostep = FALSE, save = TRUE)
fit  = gamlss(formula = Payment/Claims ~ Kilometres + Zone + Bonus + Make + Insured, Payment/Claims ~ Kilometres + Zone, family = BP(mu.link="log", sigma.link = "log"), control=conh, method=RS())

beta.est.mv = fit$mu.coefficients
nu.est.mv   = fit$sigma.coefficients

A = c(beta.est.mv, nu.est.mv)
#####################################################################
#                          RESULTS 
#####################################################################
result_est = c(beta.est.mv, nu.est.mv)

result_est
##   (Intercept)    Kilometres          Zone         Bonus          Make 
##  8.358991e+00  6.028501e-03 -1.137296e-02  7.080279e-02 -2.369009e-03 
##       Insured   (Intercept)    Kilometres          Zone 
##  7.698916e-06  1.294617e+00 -6.764708e-01 -3.025221e-01
#####################################################################
#                       PARAMETROS ESTIMADOS
#####################################################################
nEst = matrix(round(c(result_est),1), ncol = 1)
#rownames(nEst) = c("MLE")
#colnames(nEst) = c("beta0","beta1","beta2","beta3","beta4","beta5")

print.est <- function(par){
  cat("\n Beta's estimates: \n")
  print(par[1:p,])
  cat("\n Nu's estimates:\n") 
  print(par[(p+1):(p+q),])
}
print.est(nEst)
## 
##  Beta's estimates: 
## [1] 8.4 0.0 0.0 0.1 0.0 0.0
## 
##  Nu's estimates:
## [1]  1.3 -0.7 -0.3

Aplicaciones

  • \(\textbf{Análisis del riesgo-rendimiento con el modelo beta en las Cooperativas de Ahorro y Crédito (COAC's) (Andrés Felipe Sandoval Vera 2018)}\)

La implementación del Código Orgánico y Financiera a partir del año 2015, generó diferentes tipos de reacciones en sector de la economía ecuatorial sobre todo en el sector de los Seguros y de las Cooperativas de ahorro y crédito, así que proponen analizar este comportamiento en el sector cooperativo como un indicador de Beta Comparable ayuda a visualizar el panorama futuro del riesgo-rendimiento de las cooperativas analizadas frente a su mercado donde se desarrollan, mediante recopilación de datos entregados por la SEPS.

  • \(\textbf{La tasa social de descuento en la evaluación de proyectos de inversión: Aplicación para el Ecuador}\) \(\textbf{( José Gabriel Castillo y Donald Zhangallimbay 2021)}\)

En este trabajo presentan una alternativa de determinación de la tasa social de descuento basada en el modelo de estimación gamma (Weitzman, 2001), en el contexto de países en desarrollo.

Modelo Gamma

Introdución

La distribución gamma se puede utilizar con gran flexibilidad en el análisis de variables aleatorias positivas. Por lo tanto, los modelos de regresión gamma se aplican en una amplia gama de aplicaciones empíricas, como en el proceso de fijación de tasas en el marco de trabajo de carteras de seguros heterogéneas, que es la función más importante de las aseguradoras (Krishnamoorthy, 2006), y en los ingresos hospitalarios por enfermedades raras donde las series temporales son muy escasas debido a la infrecuencia de los eventos (Winklemann, 2008).

Distribución gamma

La distribución Gamma tiene la función de densidad de probabilidad:

Si \(\displaystyle{X \sim \Gamma (\alpha ,\lambda )}\) entonces su función de densidad es

\[{\displaystyle f(x)={\displaystyle\frac {\lambda (\lambda x)^{\alpha -1}e^{-\lambda x}}{\Gamma (\alpha )}}}\] para \(\displaystyle x >0\) \(\beta > 0\) es el parámetro de escala y \(\alpha > 0\) es el parámetro de forma donde: \(\Gamma (\alpha ) = \displaystyle\int_{0}^{\infty }t^{\alpha -1}e^{-t}dt\)

Modelo de regresión gamma

Sea \(Yi \sim G(\mu_i, \alpha)\), \(i = 1, . . . , n\), sea una muestra aleatoria de tamaño \(n\). En los modelos de regresión gamma con parámetro de forma constante, la estructura de regresión media se define por: \[g(\mu_i) = \eta_i = x_i'\beta\] donde \(g\) es la función de enlace, \(\beta = (\beta_0, . . ., \beta_p)'\) es el vector de media de los parámetros de regresión, \(x_i\) es el i-ésimo valor del vector de las variables explicativas y \(\eta_i\) es un predictor lineal. Algunas funciones de vínculo usuales en la regresión gamma son:

  • la función logarítmica, \(g(\mu)\) = \(log(\mu)\);
  • la función identidad, \(g(\mu) = \mu\),
  • la función inversa \(g(\mu) = 1/\mu\).

En los modelos lineales generalizados, el vínculo canónico de la media es la función inversa.

MLE

Introducción y Función Log-Verosimilitud

\[L(\alpha, \lambda) = \displaystyle\prod_{i=0}^n\left[ \displaystyle\frac{\lambda^{\alpha}e^{-\alpha x_i}(x_i)^{\alpha-1}}{\Gamma(\alpha)}\right] \] \[log(L(\alpha, \lambda)) = \log \left( \displaystyle \displaystyle\frac{\lambda^{n\alpha}e^{-\alpha \sum x_i}( \displaystyle\prod_{i=1}^n x_i)^{\alpha-1}}{[\Gamma(\alpha)]^n} \right)\] \[n\cdot\alpha log(\lambda) - \lambda\sum x_i + (\alpha - 1)log(\displaystyle\prod_{i=1}^n x_i) - log([\Gamma(\alpha)]^n) \] \[n\cdot\alpha log(\lambda) - \lambda\sum x_i + (\alpha - 1)\displaystyle\sum_{i=1}^n log(x_i) - n\cdot log(\Gamma(\alpha)) \]

MLE de \(\lambda\)

\[\displaystyle\frac{\partial l(\alpha, \lambda)}{\partial \lambda} = \displaystyle\frac{n\cdot\alpha}{\lambda} - \lambda\sum x_i = 0\] \[\displaystyle\frac{n\cdot\alpha}{\lambda} = \lambda\sum x_i \Rightarrow \lambda = \displaystyle\frac{n\cdot\alpha}{\lambda\sum x_i} \]

\[\hat{\lambda} =\displaystyle\frac{n\cdot\hat{\alpha}}{\sum x_i} = \displaystyle\frac{n\cdot\hat{\alpha}}{n \cdot \bar{x}} = \displaystyle\frac{\hat{\alpha}}{\bar{x}}\]

MLE de \(\alpha\)

\[\displaystyle\frac{\partial l(\alpha, \lambda)}{\partial \alpha} = n \cdot log(\lambda) + \displaystyle\sum_{i=1}^nlog(x_i) - \displaystyle\frac{n \cdot\Gamma'(\alpha)}{\Gamma(\alpha)} \] \[n \cdot log(\lambda) + \displaystyle\sum_{i=1}^nlog(x_i) - n \cdot DigammaFuntion(\alpha)\] \[n \cdot log(\displaystyle\frac{\hat{\alpha}}{\bar{x}}) + \displaystyle\sum_{i=1}^nlog(x_i) - n \cdot DigammaFuntion(\alpha)\] \[\hat{\alpha} = n \cdot log(\hat{\alpha}) - n \cdot log(\bar{x}) + \displaystyle\sum_{i=1}^nlog(x_i) - n \cdot DigammaFuntion(\hat{\alpha})\] \[\hat{\alpha} = n [ log(\hat{\alpha}) - log(\bar{x}) + \displaystyle\frac{1}{n}\sum_{i=1}^nlog(x_i) - DigammaFuntion(\hat{\alpha}) ]\] \[log(\hat{\alpha}) - log(\bar{x}) + \displaystyle\frac{1}{n}\sum_{i=1}^nlog(x_i) - DigammaFuntion(\hat{\alpha}) = 0\]

MLE en R

Estimaremos los \(\alpha\) y \(\lambda\).

# Establecer una seed para que el código se pueda repetir para poder reproducir resultados
set.seed(33)

# Cree datos simulados para 1000 observaciones de una variable aleatoria distribuida gamma
sampleData <- rgamma(n = 1000,shape = 2, rate =0.2)

# Sea x_bar la media muestral de los datos.
x_bar  <- mean(sampleData) # alfa/lambda es la media


# Tome el logaritmo natural de cada uno de los valores
lnSample <- log(sampleData)


# Obtener la media de los datos de muestra de registro natural
meanlnSample <- mean(lnSample)


# Escriba una función que corresponda a la derivada parcial que derivamos de la función de verosimilitud logarítmica
# con respecto a alfa.

f <- function(x){
  # En este caso de funciones, tomamos el valor x, que será el valor de alfa que estamos resolviendo.
  log(x) - digamma(x) - log(x_bar) + meanlnSample
}

# Use la función uniroot para encontrar dónde la derivada es igual a cero para que podamos encontrar
# la raíz (valor x o sombrero alfa) en la que ocurre.
a <- uniroot(f,lower = 2, upper = 2.2)

# Establecer alpha_hat igual a la raíz de la función f en el intervalo (2.0,2.2)
alpha_hat <- a$root
alpha_hat
## [1] 2.080856
# Recuerde que se encontró que la estimación de probabilidad máxima de Lambda era alpha_hat sobre x_bar.
lambda_hat <- alpha_hat/x_bar

lambda_hat
## [1] 0.2096357

Aplicaciones

  • \(\textbf{Modelación de la distribución gamma en matlab para aplicaciones de radar ( José Raul Machado-Fernandez 2016)}\)

El clutter es una señal de interferencia que aparece en los sistemas de radar. El autor modeló la distribución Gamma, íntimamente relacionado con clutter, en MATLAB. La implementación permite un acceso fácil a la manipulación de las funciones de densidad de probabilidad y de distribución, la generación de muestras, el cálculo de los momentos, la aplicación de algoritmos de bondad de ajuste y la estimación de parámetros a partir de muestras con distribución gamma.

  • \(\textbf{La tasa social de descuento en la evaluación de proyectos de inversión: Aplicación para el Ecuador}\) \(\textbf{( José Gabriel Castillo y Donald Zhangallimbay 2021)}\)

En este trabajo presentan una alternativa de determinación de la tasa social de descuento basada en el modelo de estimación gamma (Weitzman, 2001), en el contexto de países en desarrollo.

  • \(\textbf{Reclamaciones de seguros}\)

La distribución gamma se ha utilizado para modelar el tamaño de las reclamaciones de seguros\(^{[1]}\) y las precipitaciones\(^{[2]}\). Esto significa que las reclamaciones de seguros agregadas y la cantidad de lluvia acumulada en un embalse se modelan mediante un proceso gamma , al igual que la distribución exponencial genera un proceso de Poisson.

  • \(\textbf{Comunicación Inalámbrica}\)

La distribución gamma se utiliza para modelar el desvanecimiento de la potencia de la señal por múltiples rutas.

  • \(\textbf{Oncología}\)

La distribución por edades de la incidencia de cáncer a menudo sigue la distribución gamma, mientras que los parámetros de forma y escala predicen, respectivamente, el número de eventos impulsores y el intervalo de tiempo entre ellos\(^{[3]}\).

  • \(\textbf{Neurociencia}\)

la distribución gamma se usa a menudo para describir la distribución de intervalos entre picos\(^{[4][5]}\).

  • \(\textbf{Expresión génica bacteriana}\)

El número de copias de una proteína expresada constitutivamente a menudo sigue la distribución gamma, donde el parámetro de escala y forma son, respectivamente, el número medio de explosiones por ciclo celular y el número medio de moléculas de proteína producidas por un solo ARNm durante su vida\(^{[6]}\).

  • \(\textbf{En genómica}\)

La distribución gamma se aplicó en el paso de llamada de pico (es decir, en el reconocimiento de la señal) en el análisis de datos de ChIP-chip\(^{[7]}\) y ChIP-seq\(^{[8]}\).

  • \(\textbf{Estadística bayesiana}\)

La distribución gamma se usa ampliamente como un previo conjugado en la estadística bayesiana. Es el conjugado previo para la precisión (es decir, inversa de la varianza) de una distribución normal. También es el conjugado previo de la distribución exponencial.

Notas

  1. p. 43, Philip J. Boland, Métodos estadísticos y probabilísticos en ciencia actuarial, Chapman & Hall CRC 2007

  2. Aksoy, H. (2000) “Uso de distribución gamma en análisis hidrológico” , Turk J. Engin Environ Sci , 24, 419 - 428.

  3. Belikov, Aleksey V. (22 de septiembre de 2017). “El número de eventos cancerígenos clave se puede predecir a partir de la incidencia del cáncer” . Informes científicos. 7 (1): 12170. doi : 10.1038 / s41598-017-12448-7. PMC 5610194 . PMID 28939880 .

  4. JG Robson y JB Troy, “Naturaleza de la descarga mantenida de las células ganglionares retinianas Q, X e Y del gato”, J. Opt. Soc. Soy. A 4, 2301-2307 (1987)

  5. MCM Wright, IM Winter, JJ Forster, S. Bleeck “La respuesta a las ráfagas de tono de la mejor frecuencia en el núcleo coclear ventral se rige por las estadísticas ordenadas del intervalo entre picos”, Hearing Research 317 (2014)

  6. N. Friedman, L. Cai y XS Xie (2006) “Vinculación de la dinámica estocástica con la distribución de la población: un marco analítico de expresión génica”, Phys. Rev. Lett. 97, 168302.

  7. DJ Reiss, MT Facciotti y NS Baliga (2008) “Deconvolución basada en modelos de la unión del ADN en todo el genoma” , Bioinformática , 24, 396–403.

  8. MA Mendoza-Parra, M Nowicka, W Van Gool, H Gronemeyer (2013) “Caracterización de patrones de unión de ChIP-seq por deconvolución de forma de pico basada en modelo” , BMC Genomics , 14: 834

Modelos

library(tidyverse)

csv_data <- read_csv('midata.xls')

head(csv_data)
## # A tibble: 6 × 8
##    ...1 Kilometres  Zone Bonus  Make Insured Claims Payment
##   <dbl>      <dbl> <dbl> <dbl> <dbl>   <dbl>  <dbl>   <dbl>
## 1     1          1     1     1     1   455.     108  392491
## 2     2          1     1     1     2    69.2     19   46221
## 3     3          1     1     1     3    72.9     13   15694
## 4     4          1     1     1     4  1292.     124  422201
## 5     5          1     1     1     5   191.      40  119373
## 6     6          1     1     1     6   478.      57  170913
# Datos Corregidos
csv_Data <- csv_data[2:8]

Kilometres <-  csv_Data$Kilometres
Zone <- csv_Data$Zone
Bonus <- csv_Data$Bonus
Make <- csv_Data$Make
Insured <- csv_Data$Insured
Claims <- csv_Data$Claims
Payment <- csv_Data$Payment

head(csv_Data)
## # A tibble: 6 × 7
##   Kilometres  Zone Bonus  Make Insured Claims Payment
##        <dbl> <dbl> <dbl> <dbl>   <dbl>  <dbl>   <dbl>
## 1          1     1     1     1   455.     108  392491
## 2          1     1     1     2    69.2     19   46221
## 3          1     1     1     3    72.9     13   15694
## 4          1     1     1     4  1292.     124  422201
## 5          1     1     1     5   191.      40  119373
## 6          1     1     1     6   478.      57  170913

Modelo M1

Este modelo tomara todas las variables de la base de datos.

m1 = glm(formula = Payment/Claims ~ Kilometres + Zone + Bonus + Make + Insured , family =Gamma(link="inverse"),
         data = csv_Data)
summary(m1)
## 
## Call:
## glm(formula = Payment/Claims ~ Kilometres + Zone + Bonus + Make + 
##     Insured, family = Gamma(link = "inverse"), data = csv_Data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5885  -0.6075  -0.1533   0.1478   2.6877  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.387e-04  1.581e-05  15.102  < 2e-16 ***
## Kilometres  -5.210e-06  2.772e-06  -1.880  0.06033 .  
## Zone        -5.426e-06  2.068e-06  -2.624  0.00876 ** 
## Bonus       -1.756e-06  1.911e-06  -0.919  0.35836    
## Make        -8.071e-07  1.459e-06  -0.553  0.58009    
## Insured      1.367e-10  6.766e-10   0.202  0.83996    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.7162004)
## 
##     Null deviance: 899.37  on 1796  degrees of freedom
## Residual deviance: 890.94  on 1791  degrees of freedom
##   (385 observations deleted due to missingness)
## AIC: 33857
## 
## Number of Fisher Scoring iterations: 6

Modelo M2

Criterio de Selección

El método utilizado es el método forward. Para aplicar este método se debe crear un modelo vacío (empty.model) del cual iniciará el proceso. Es necesario definir un punto final de búsqueda, ese punto es una fórmula que en este caso seran todas las variables.

Código en R

Se usa la función stepAIC (Extraída de la librería MASS) y se elije trace=TRUE para que SI se muestren los detalles del proceso de selección.

Proceso de selección:

empty.model <-  glm(formula = Payment/Claims ~ 1,data = csv_Data,family = Gamma(link="inverse") )

m2 <- stepAIC(empty.model, trace=TRUE, direction="forward", scope = Payment/Claims ~ Kilometres + Zone + Bonus + Make + Insured)
## Start:  AIC=33865.11
## Payment/Claims ~ 1
## 
##              Df Deviance   AIC
## + Zone        1   894.48 33861
## + Kilometres  1   896.94 33864
## <none>            899.37 33865
## + Bonus       1   898.57 33866
## + Make        1   899.09 33867
## + Insured     1   899.33 33867
## 
## Step:  AIC=33856.5
## Payment/Claims ~ Zone
## 
##              Df Deviance   AIC
## + Kilometres  1   891.74 33855
## <none>            894.48 33857
## + Bonus       1   893.83 33858
## + Make        1   894.27 33858
## + Insured     1   894.46 33858
## 
## Step:  AIC=33852.54
## Payment/Claims ~ Zone + Kilometres
## 
##           Df Deviance   AIC
## <none>         891.74 33853
## + Bonus    1   891.17 33854
## + Make     1   891.55 33854
## + Insured  1   891.74 33855

Para obtener un resumen del proceso se usa:

m2$anova
## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## Payment/Claims ~ 1
## 
## Final Model:
## Payment/Claims ~ Zone + Kilometres
## 
## 
##           Step Df Deviance Resid. Df Resid. Dev      AIC
## 1                               1796   899.3748 33865.11
## 2       + Zone  1 4.896303      1795   894.4785 33856.50
## 3 + Kilometres  1 2.739782      1794   891.7387 33852.54

Para ver la tabla de resultados del modelo m2.

summary(m2)
## 
## Call:
## glm(formula = Payment/Claims ~ Zone + Kilometres, family = Gamma(link = "inverse"), 
##     data = csv_Data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5936  -0.6107  -0.1564   0.1458   2.6953  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.286e-04  1.217e-05  18.782  < 2e-16 ***
## Zone        -5.574e-06  2.055e-06  -2.713  0.00673 ** 
## Kilometres  -5.384e-06  2.738e-06  -1.967  0.04939 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.71051)
## 
##     Null deviance: 899.37  on 1796  degrees of freedom
## Residual deviance: 891.74  on 1794  degrees of freedom
##   (385 observations deleted due to missingness)
## AIC: 33853
## 
## Number of Fisher Scoring iterations: 6

Modelo M3

Criterio de Selección

Para aplicar este método se debe crear un modelo vacío del cual iniciará el proceso. Es necesario definir un punto final de búsqueda, ese punto es una formula que en este caso llamaremos horizonte. A continuación el código.

Código en R

Se usa la función stepAIC (Extraída de la librería MASS) y se elije trace=TRUE para que SI se muestren los detalles del proceso de selección.

Proceso de selección:

empty.model <- glm(formula = Payment/Claims ~1,data = csv_Data,family = Gamma(link="inverse") )


m3 <- stepAIC(empty.model, trace=TRUE, direction="both", scope = Payment/Claims ~ Kilometres + Zone + Bonus + Make + Insured)
## Start:  AIC=33865.11
## Payment/Claims ~ 1
## 
##              Df Deviance   AIC
## + Zone        1   894.48 33861
## + Kilometres  1   896.94 33864
## <none>            899.37 33865
## + Bonus       1   898.57 33866
## + Make        1   899.09 33867
## + Insured     1   899.33 33867
## 
## Step:  AIC=33856.5
## Payment/Claims ~ Zone
## 
##              Df Deviance   AIC
## + Kilometres  1   891.74 33855
## <none>            894.48 33857
## + Bonus       1   893.83 33858
## + Make        1   894.27 33858
## + Insured     1   894.46 33858
## - Zone        1   899.37 33861
## 
## Step:  AIC=33852.54
## Payment/Claims ~ Zone + Kilometres
## 
##              Df Deviance   AIC
## <none>            891.74 33853
## + Bonus       1   891.17 33854
## + Make        1   891.55 33854
## - Kilometres  1   894.48 33854
## + Insured     1   891.74 33855
## - Zone        1   896.94 33858

Para obtener un resumen del proceso se usa:

m3$anova
## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## Payment/Claims ~ 1
## 
## Final Model:
## Payment/Claims ~ Zone + Kilometres
## 
## 
##           Step Df Deviance Resid. Df Resid. Dev      AIC
## 1                               1796   899.3748 33865.11
## 2       + Zone  1 4.896303      1795   894.4785 33856.50
## 3 + Kilometres  1 2.739782      1794   891.7387 33852.54

Para ver la tabla de resultados del modelo m3.

summary(m3)
## 
## Call:
## glm(formula = Payment/Claims ~ Zone + Kilometres, family = Gamma(link = "inverse"), 
##     data = csv_Data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5936  -0.6107  -0.1564   0.1458   2.6953  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.286e-04  1.217e-05  18.782  < 2e-16 ***
## Zone        -5.574e-06  2.055e-06  -2.713  0.00673 ** 
## Kilometres  -5.384e-06  2.738e-06  -1.967  0.04939 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.71051)
## 
##     Null deviance: 899.37  on 1796  degrees of freedom
## Residual deviance: 891.74  on 1794  degrees of freedom
##   (385 observations deleted due to missingness)
## AIC: 33853
## 
## Number of Fisher Scoring iterations: 6

Modelo M4

conh = gamlss.control(trace = FALSE, autostep = FALSE, save = TRUE)
m4 <- gamlss(Payment/Claims ~ Kilometres + Zone + Bonus + Make + Insured, family = BP(mu.link="log", sigma.link = "log"), control=conh, method=RS())

summary(m4)
## Warning in summary.gamlss(m4): summary: vcov has failed, option qr is used instead
## ******************************************************************
## Family:  c("BP", "Beta Prime") 
## 
## Call:  gamlss(formula = Payment/Claims ~ Kilometres + Zone + Bonus +  
##     Make + Insured, family = BP(mu.link = "log", sigma.link = "log"),  
##     method = RS(), control = conh) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.692e+00  6.717e-02 129.392  < 2e-16 ***
## Kilometres  -4.867e-02  1.207e-02  -4.033 5.74e-05 ***
## Zone        -3.754e-02  9.024e-03  -4.160 3.33e-05 ***
## Bonus        7.348e-02  8.297e-03   8.856  < 2e-16 ***
## Make        -5.003e-04  6.362e-03  -0.079   0.9373    
## Insured      7.703e-06  2.794e-06   2.757   0.0059 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -4.201      1.916  -2.192   0.0285 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1797 
## Degrees of Freedom for the fit:  7
##       Residual Deg. of Freedom:  1790 
##                       at cycle:  20 
##  
## Global Deviance:     33992.78 
##             AIC:     34006.78 
##             SBC:     34045.24 
## ******************************************************************

Comparaciones de los Modelos

Comparación por Residuos

par(mfrow=c(1, 4))
plot(m1, main="M1", pch=19, cex=1, which=1)
plot(m2, main="M2 Forward", pch=19, cex=1, which=1)
plot(m3, main="M3 Both", pch=19, cex=1, which=1)
plot(m4, main="M4", pch=19, cex=1, which=1)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  0.01601973 
##                        variance   =  0.993483 
##                coef. of skewness  =  -1.648564 
##                coef. of kurtosis  =  13.8506 
## Filliben correlation coefficient  =  0.957701 
## ******************************************************************

Conclusiones por Residuos

  • Notamos que para los modelos gamma los presentaban grandes similitudes por lo tanto no es claro concluir que modelo gamma es mejor.

  • El modelo con familia BetaPrime es explicativo pero presenta mejor ajuste los tres modelos gamma.