Discrete power-Ailamujia distribution

Source: R/DsPA.R

The function DsPA() defines the discrete power-Ailamujia distribution a two parameter distribution, for a gamlss.family object to be used in GAMLSS fitting using the function gamlss().

DsPA(mu.link = "log", sigma.link = "logit")

Arguments

defines the mu.link, with "log" link as the default for the mu parameter.

defines the sigma.link, with "logit" link as the default for the sigma parameter. Other links are "probit" and "cloglog"'(complementary log-log).

Value

Returns a gamlss.family object which can be used to fit a discrete power-Ailamujia distribution in the gamlss() function.

Details

The discrete power-Ailamujia distribution with parameters \(\mu\) and \(\sigma\) has a support 0, 1, 2, ... and density given by

\(f(x | \mu, \sigma) = (\sigma^x)^\mu (1-x^\mu \log(\lambda)) - (\sigma^{(x+1)})^\mu (1-(x+1)^\mu \log(\lambda))\)

Note: in this implementation we changed the original parameters \(\beta\) and \(\lambda\) for \(\mu\) and \(\sigma\) respectively, we did it to implement this distribution within gamlss framework.

References

Alghamdi, A. S., Ahsan-ul-Haq, M., Babar, A., Aljohani, H. M., Afify, A. Z., & Cell, Q. E. (2022). The discrete power-Ailamujia distribution: properties, inference, and applications. AIMS Math, 7(5), 8344-8360.

See also

dDsPA.

Author

Maria Camila Mena Romana, mamenar@unal.edu.co

Examples

# Example 1
# Generating some random values with
# known mu and sigma
y <- rDsPA(n=100, mu=1.2, sigma=0.5)

# Fitting the model
library(gamlss)
mod1 <- gamlss(y~1, family=DsPA,
               control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 340.991 

# Extracting the fitted values for mu and sigma
# using the inverse link function
inv_logit <- function(x) 1/(1 + exp(-x))

exp(coef(mod1, what="mu"))
#> (Intercept) 
#>    1.220429 
inv_logit(coef(mod1, what="sigma"))
#> (Intercept) 
#>   0.5264461 

# Example 2
# Generating random values under some model

# A function to simulate a data set with Y ~ GGEO
gendat <- function(n) {
  x1 <- runif(n)
  x2 <- runif(n)
  x3 <- runif(n)
  x4 <- runif(n)
  mu    <- exp(1 + 1.2*x1 + 0.2*x2)
  sigma <- inv_logit(2 + 1.5*x3 + 1.5*x4)
  y <- rDsPA(n=n, mu=mu, sigma=sigma)
  data.frame(y=y, x1=x1, x2=x2, x3=x3, x4=x4)
}

set.seed(123)
datos <- gendat(n=100)

mod2 <- gamlss(y~x1+x2, sigma.fo=~x3+x4, family=DsPA, data=datos,
               control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 167.4707 
#> GAMLSS-RS iteration 2: Global Deviance = 158.9519 
#> GAMLSS-RS iteration 3: Global Deviance = 152.8383 
#> GAMLSS-RS iteration 4: Global Deviance = 148.164 
#> GAMLSS-RS iteration 5: Global Deviance = 144.4253 
#> GAMLSS-RS iteration 6: Global Deviance = 141.393 
#> GAMLSS-RS iteration 7: Global Deviance = 138.9571 
#> GAMLSS-RS iteration 8: Global Deviance = 136.9419 
#> GAMLSS-RS iteration 9: Global Deviance = 135.2594 
#> GAMLSS-RS iteration 10: Global Deviance = 133.8517 
#> GAMLSS-RS iteration 11: Global Deviance = 132.6168 
#> GAMLSS-RS iteration 12: Global Deviance = 131.5779 
#> GAMLSS-RS iteration 13: Global Deviance = 130.656 
#> GAMLSS-RS iteration 14: Global Deviance = 129.844 
#> GAMLSS-RS iteration 15: Global Deviance = 129.1531 
#> GAMLSS-RS iteration 16: Global Deviance = 128.5258 
#> GAMLSS-RS iteration 17: Global Deviance = 127.9661 
#> GAMLSS-RS iteration 18: Global Deviance = 127.4814 
#> GAMLSS-RS iteration 19: Global Deviance = 127.0238 
#> GAMLSS-RS iteration 20: Global Deviance = 126.6203 
#> GAMLSS-RS iteration 21: Global Deviance = 126.2604 
#> GAMLSS-RS iteration 22: Global Deviance = 125.9382 
#> GAMLSS-RS iteration 23: Global Deviance = 125.6367 
#> GAMLSS-RS iteration 24: Global Deviance = 125.3682 
#> GAMLSS-RS iteration 25: Global Deviance = 125.1274 
#> GAMLSS-RS iteration 26: Global Deviance = 124.9254 
#> GAMLSS-RS iteration 27: Global Deviance = 124.7293 
#> GAMLSS-RS iteration 28: Global Deviance = 124.567 
#> GAMLSS-RS iteration 29: Global Deviance = 124.4201 
#> GAMLSS-RS iteration 30: Global Deviance = 124.2733 
#> GAMLSS-RS iteration 31: Global Deviance = 124.1475 
#> GAMLSS-RS iteration 32: Global Deviance = 124.0343 
#> GAMLSS-RS iteration 33: Global Deviance = 123.9309 
#> GAMLSS-RS iteration 34: Global Deviance = 123.8347 
#> GAMLSS-RS iteration 35: Global Deviance = 123.747 
#> GAMLSS-RS iteration 36: Global Deviance = 123.6644 
#> GAMLSS-RS iteration 37: Global Deviance = 123.5897 
#> GAMLSS-RS iteration 38: Global Deviance = 123.5298 
#> GAMLSS-RS iteration 39: Global Deviance = 123.4714 
#> GAMLSS-RS iteration 40: Global Deviance = 123.4117 
#> GAMLSS-RS iteration 41: Global Deviance = 123.3276 
#> GAMLSS-RS iteration 42: Global Deviance = 123.2816 
#> GAMLSS-RS iteration 43: Global Deviance = 123.2117 
#> GAMLSS-RS iteration 44: Global Deviance = 123.1678 
#> GAMLSS-RS iteration 45: Global Deviance = 123.111 
#> GAMLSS-RS iteration 46: Global Deviance = 123.0746 
#> GAMLSS-RS iteration 47: Global Deviance = 123.0289 
#> GAMLSS-RS iteration 48: Global Deviance = 122.9986 
#> GAMLSS-RS iteration 49: Global Deviance = 122.9843 
#> GAMLSS-RS iteration 50: Global Deviance = 122.9627 
#> GAMLSS-RS iteration 51: Global Deviance = 122.9393 
#> GAMLSS-RS iteration 52: Global Deviance = 122.9263 
#> GAMLSS-RS iteration 53: Global Deviance = 122.9165 
#> GAMLSS-RS iteration 54: Global Deviance = 122.8971 
#> GAMLSS-RS iteration 55: Global Deviance = 122.8806 
#> GAMLSS-RS iteration 56: Global Deviance = 122.8695 
#> GAMLSS-RS iteration 57: Global Deviance = 122.8605 
#> GAMLSS-RS iteration 58: Global Deviance = 122.8528 
#> GAMLSS-RS iteration 59: Global Deviance = 122.8461 
#> GAMLSS-RS iteration 60: Global Deviance = 122.8399 
#> GAMLSS-RS iteration 61: Global Deviance = 122.8341 
#> GAMLSS-RS iteration 62: Global Deviance = 122.8284 
#> GAMLSS-RS iteration 63: Global Deviance = 122.8229 
#> GAMLSS-RS iteration 64: Global Deviance = 122.8174 
#> GAMLSS-RS iteration 65: Global Deviance = 122.812 
#> GAMLSS-RS iteration 66: Global Deviance = 122.8067 
#> GAMLSS-RS iteration 67: Global Deviance = 122.8014 
#> GAMLSS-RS iteration 68: Global Deviance = 122.7961 
#> GAMLSS-RS iteration 69: Global Deviance = 122.7909 
#> GAMLSS-RS iteration 70: Global Deviance = 122.7847 
#> GAMLSS-RS iteration 71: Global Deviance = 122.7789 
#> GAMLSS-RS iteration 72: Global Deviance = 122.7731 
#> GAMLSS-RS iteration 73: Global Deviance = 122.7675 
#> GAMLSS-RS iteration 74: Global Deviance = 122.7619 
#> GAMLSS-RS iteration 75: Global Deviance = 122.7565 
#> GAMLSS-RS iteration 76: Global Deviance = 122.7511 
#> GAMLSS-RS iteration 77: Global Deviance = 122.7459 
#> GAMLSS-RS iteration 78: Global Deviance = 122.7408 
#> GAMLSS-RS iteration 79: Global Deviance = 122.7358 
#> GAMLSS-RS iteration 80: Global Deviance = 122.7309 
#> GAMLSS-RS iteration 81: Global Deviance = 122.7261 
#> GAMLSS-RS iteration 82: Global Deviance = 122.7215 
#> GAMLSS-RS iteration 83: Global Deviance = 122.717 
#> GAMLSS-RS iteration 84: Global Deviance = 122.7126 
#> GAMLSS-RS iteration 85: Global Deviance = 122.7083 
#> GAMLSS-RS iteration 86: Global Deviance = 122.7042 
#> GAMLSS-RS iteration 87: Global Deviance = 122.7002 
#> GAMLSS-RS iteration 88: Global Deviance = 122.6963 
#> GAMLSS-RS iteration 89: Global Deviance = 122.6925 
#> GAMLSS-RS iteration 90: Global Deviance = 122.6889 
#> GAMLSS-RS iteration 91: Global Deviance = 122.6853 
#> GAMLSS-RS iteration 92: Global Deviance = 122.6819 
#> GAMLSS-RS iteration 93: Global Deviance = 122.6786 
#> GAMLSS-RS iteration 94: Global Deviance = 122.6755 
#> GAMLSS-RS iteration 95: Global Deviance = 122.6724 
#> GAMLSS-RS iteration 96: Global Deviance = 122.6695 
#> GAMLSS-RS iteration 97: Global Deviance = 122.6666 
#> GAMLSS-RS iteration 98: Global Deviance = 122.6639 
#> GAMLSS-RS iteration 99: Global Deviance = 122.6612 
#> GAMLSS-RS iteration 100: Global Deviance = 122.6587 
#> GAMLSS-RS iteration 101: Global Deviance = 122.6563 
#> GAMLSS-RS iteration 102: Global Deviance = 122.6539 
#> GAMLSS-RS iteration 103: Global Deviance = 122.6517 
#> GAMLSS-RS iteration 104: Global Deviance = 122.6495 
#> GAMLSS-RS iteration 105: Global Deviance = 122.6474 
#> GAMLSS-RS iteration 106: Global Deviance = 122.6455 
#> GAMLSS-RS iteration 107: Global Deviance = 122.6435 
#> GAMLSS-RS iteration 108: Global Deviance = 122.6417 
#> GAMLSS-RS iteration 109: Global Deviance = 122.64 
#> GAMLSS-RS iteration 110: Global Deviance = 122.6383 
#> GAMLSS-RS iteration 111: Global Deviance = 122.6367 
#> GAMLSS-RS iteration 112: Global Deviance = 122.6351 
#> GAMLSS-RS iteration 113: Global Deviance = 122.6337 
#> GAMLSS-RS iteration 114: Global Deviance = 122.6323 
#> GAMLSS-RS iteration 115: Global Deviance = 122.6309 
#> GAMLSS-RS iteration 116: Global Deviance = 122.6296 
#> GAMLSS-RS iteration 117: Global Deviance = 122.6267 
#> GAMLSS-RS iteration 118: Global Deviance = 122.6248 
#> GAMLSS-RS iteration 119: Global Deviance = 122.6234 
#> GAMLSS-RS iteration 120: Global Deviance = 122.6214 
#> GAMLSS-RS iteration 121: Global Deviance = 122.6203 
#> GAMLSS-RS iteration 122: Global Deviance = 122.6186 
#> GAMLSS-RS iteration 123: Global Deviance = 122.6173 
#> GAMLSS-RS iteration 124: Global Deviance = 122.616 
#> GAMLSS-RS iteration 125: Global Deviance = 122.6149 
#> GAMLSS-RS iteration 126: Global Deviance = 122.6139 
#> GAMLSS-RS iteration 127: Global Deviance = 122.6129 

summary(mod2)
#> Warning: summary: vcov has failed, option qr is used instead
#> ******************************************************************
#> Family:  c("DsPA", "discrete power-Ailamujia") 
#> 
#> Call:  gamlss(formula = y ~ x1 + x2, sigma.formula = ~x3 + x4, family = DsPA,  
#>     data = datos, control = gamlss.control(n.cyc = 500, trace = TRUE)) 
#> 
#> Fitting method: RS() 
#> 
#> ------------------------------------------------------------------
#> Mu link function:  log
#> Mu Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   0.9618     0.0880  10.929   <2e-16 ***
#> x1            1.2432     0.1244   9.993   <2e-16 ***
#> x2            0.1302     0.1296   1.004    0.318    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> Sigma link function:  logit
#> Sigma Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   1.4613     0.3826   3.820 0.000236 ***
#> x3            1.4847     0.4191   3.542 0.000611 ***
#> x4            2.1236     0.4093   5.189 1.16e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> No. of observations in the fit:  100 
#> Degrees of Freedom for the fit:  6
#>       Residual Deg. of Freedom:  94 
#>                       at cycle:  127 
#>  
#> Global Deviance:     122.6129 
#>             AIC:     134.6129 
#>             SBC:     150.244 
#> ******************************************************************

# Example 3
# failure times for a sample of 15 electronic components in an acceleration life test
# Taken from
# Alghamdi et. al (202) page 8354.

y <- c(1.0, 5.0, 6.0, 11.0, 12.0, 19.0, 20.0, 22.0,
       23.0, 31.0, 37.0, 46.0, 54.0, 60.0, 66.0)

mod3 <- gamlss(y~1, family=DsPA,
               control=gamlss.control(n.cyc=500, trace=FALSE))

# Extracting the fitted values for mu and sigma
# using the inverse link function
inv_logit <- function(x) 1/(1 + exp(-x))

exp(coef(mod3, what="mu"))
#> (Intercept) 
#>   0.8590179 
inv_logit(coef(mod3, what="sigma"))
#> (Intercept) 
#>   0.8886351 

# Example 4
# number of fires in Greece from July 1, 1998 to August 31, 1998.
# Taken from
# Alghamdi et. al (202) page 8354.

y <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3,
       3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5,
       5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6,
       6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7,
       8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9,
       9, 9, 9, 9, 10, 10, 10, 11, 11,
       11, 11, 12, 12, 12, 12, 12, 12,
       15, 15, 15, 15, 16, 20, 43)

mod4 <- gamlss(y~1, family=DsPA,
               control=gamlss.control(n.cyc=500, trace=FALSE))

# Extracting the fitted values for mu and sigma
# using the inverse link function
inv_logit <- function(x) 1/(1 + exp(-x))

exp(coef(mod4, what="mu"))
#> (Intercept) 
#>    0.770927 
inv_logit(coef(mod4, what="sigma"))
#> (Intercept) 
#>    0.581188