Discrete power-Ailamujia distribution — DsPA • DiscreteDists

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

mu.link: defines the mu.link, with "log" link as the default for the mu parameter.
sigma.link: 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

Author

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

Examples

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

# Fitting the model
library(gamlss)
mod1 <- 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(mod1, what="mu"))
#> (Intercept) 
#>    1.194329 
inv_logit(coef(mod1, what="sigma"))
#> (Intercept) 
#>   0.4975115 

# 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(16494786)
datos <- gendat(n=500)

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 = 806.0802 
#> GAMLSS-RS iteration 2: Global Deviance = 767.9744 
#> GAMLSS-RS iteration 3: Global Deviance = 739.4638 
#> GAMLSS-RS iteration 4: Global Deviance = 717.5791 
#> GAMLSS-RS iteration 5: Global Deviance = 700.2108 
#> GAMLSS-RS iteration 6: Global Deviance = 686.31 
#> GAMLSS-RS iteration 7: Global Deviance = 674.8226 
#> GAMLSS-RS iteration 8: Global Deviance = 665.1847 
#> GAMLSS-RS iteration 9: Global Deviance = 657.097 
#> GAMLSS-RS iteration 10: Global Deviance = 650.1513 
#> GAMLSS-RS iteration 11: Global Deviance = 644.2052 
#> GAMLSS-RS iteration 12: Global Deviance = 639.0427 
#> GAMLSS-RS iteration 13: Global Deviance = 634.5428 
#> GAMLSS-RS iteration 14: Global Deviance = 630.6079 
#> GAMLSS-RS iteration 15: Global Deviance = 627.2147 
#> GAMLSS-RS iteration 16: Global Deviance = 624.1372 
#> GAMLSS-RS iteration 17: Global Deviance = 621.4152 
#> GAMLSS-RS iteration 18: Global Deviance = 619.0256 
#> GAMLSS-RS iteration 19: Global Deviance = 616.7777 
#> GAMLSS-RS iteration 20: Global Deviance = 614.8309 
#> GAMLSS-RS iteration 21: Global Deviance = 613.1152 
#> GAMLSS-RS iteration 22: Global Deviance = 611.5554 
#> GAMLSS-RS iteration 23: Global Deviance = 610.1376 
#> GAMLSS-RS iteration 24: Global Deviance = 608.8251 
#> GAMLSS-RS iteration 25: Global Deviance = 607.6118 
#> GAMLSS-RS iteration 26: Global Deviance = 606.5167 
#> GAMLSS-RS iteration 27: Global Deviance = 605.547 
#> GAMLSS-RS iteration 28: Global Deviance = 604.6625 
#> GAMLSS-RS iteration 29: Global Deviance = 603.8544 
#> GAMLSS-RS iteration 30: Global Deviance = 603.1155 
#> GAMLSS-RS iteration 31: Global Deviance = 602.4547 
#> GAMLSS-RS iteration 32: Global Deviance = 601.8468 
#> GAMLSS-RS iteration 33: Global Deviance = 601.287 
#> GAMLSS-RS iteration 34: Global Deviance = 600.7715 
#> GAMLSS-RS iteration 35: Global Deviance = 600.2963 
#> GAMLSS-RS iteration 36: Global Deviance = 599.8583 
#> GAMLSS-RS iteration 37: Global Deviance = 599.4468 
#> GAMLSS-RS iteration 38: Global Deviance = 599.0683 
#> GAMLSS-RS iteration 39: Global Deviance = 598.7198 
#> GAMLSS-RS iteration 40: Global Deviance = 598.3986 
#> GAMLSS-RS iteration 41: Global Deviance = 598.1023 
#> GAMLSS-RS iteration 42: Global Deviance = 597.8288 
#> GAMLSS-RS iteration 43: Global Deviance = 597.5687 
#> GAMLSS-RS iteration 44: Global Deviance = 597.3292 
#> GAMLSS-RS iteration 45: Global Deviance = 597.1086 
#> GAMLSS-RS iteration 46: Global Deviance = 596.9052 
#> GAMLSS-RS iteration 47: Global Deviance = 596.7174 
#> GAMLSS-RS iteration 48: Global Deviance = 596.544 
#> GAMLSS-RS iteration 49: Global Deviance = 596.3837 
#> GAMLSS-RS iteration 50: Global Deviance = 596.2355 
#> GAMLSS-RS iteration 51: Global Deviance = 596.0984 
#> GAMLSS-RS iteration 52: Global Deviance = 595.9715 
#> GAMLSS-RS iteration 53: Global Deviance = 595.8568 
#> GAMLSS-RS iteration 54: Global Deviance = 595.7502 
#> GAMLSS-RS iteration 55: Global Deviance = 595.6511 
#> GAMLSS-RS iteration 56: Global Deviance = 595.5589 
#> GAMLSS-RS iteration 57: Global Deviance = 595.4731 
#> GAMLSS-RS iteration 58: Global Deviance = 595.3933 
#> GAMLSS-RS iteration 59: Global Deviance = 595.319 
#> GAMLSS-RS iteration 60: Global Deviance = 595.2498 
#> GAMLSS-RS iteration 61: Global Deviance = 595.1853 
#> GAMLSS-RS iteration 62: Global Deviance = 595.1253 
#> GAMLSS-RS iteration 63: Global Deviance = 595.0694 
#> GAMLSS-RS iteration 64: Global Deviance = 595.0173 
#> GAMLSS-RS iteration 65: Global Deviance = 594.9712 
#> GAMLSS-RS iteration 66: Global Deviance = 594.9281 
#> GAMLSS-RS iteration 67: Global Deviance = 594.8878 
#> GAMLSS-RS iteration 68: Global Deviance = 594.85 
#> GAMLSS-RS iteration 69: Global Deviance = 594.8146 
#> GAMLSS-RS iteration 70: Global Deviance = 594.7814 
#> GAMLSS-RS iteration 71: Global Deviance = 594.7503 
#> GAMLSS-RS iteration 72: Global Deviance = 594.7212 
#> GAMLSS-RS iteration 73: Global Deviance = 594.6938 
#> GAMLSS-RS iteration 74: Global Deviance = 594.6682 
#> GAMLSS-RS iteration 75: Global Deviance = 594.6442 
#> GAMLSS-RS iteration 76: Global Deviance = 594.6217 
#> GAMLSS-RS iteration 77: Global Deviance = 594.6006 
#> GAMLSS-RS iteration 78: Global Deviance = 594.5809 
#> GAMLSS-RS iteration 79: Global Deviance = 594.5623 
#> GAMLSS-RS iteration 80: Global Deviance = 594.5449 
#> GAMLSS-RS iteration 81: Global Deviance = 594.5286 
#> GAMLSS-RS iteration 82: Global Deviance = 594.5133 
#> GAMLSS-RS iteration 83: Global Deviance = 594.4989 
#> GAMLSS-RS iteration 84: Global Deviance = 594.4854 
#> GAMLSS-RS iteration 85: Global Deviance = 594.4728 
#> GAMLSS-RS iteration 86: Global Deviance = 594.4609 
#> GAMLSS-RS iteration 87: Global Deviance = 594.4497 
#> GAMLSS-RS iteration 88: Global Deviance = 594.4393 
#> GAMLSS-RS iteration 89: Global Deviance = 594.4295 
#> GAMLSS-RS iteration 90: Global Deviance = 594.4202 
#> GAMLSS-RS iteration 91: Global Deviance = 594.4116 
#> GAMLSS-RS iteration 92: Global Deviance = 594.4034 
#> GAMLSS-RS iteration 93: Global Deviance = 594.3958 
#> GAMLSS-RS iteration 94: Global Deviance = 594.3886 
#> GAMLSS-RS iteration 95: Global Deviance = 594.3819 
#> GAMLSS-RS iteration 96: Global Deviance = 594.3756 
#> GAMLSS-RS iteration 97: Global Deviance = 594.3696 
#> GAMLSS-RS iteration 98: Global Deviance = 594.364 
#> GAMLSS-RS iteration 99: Global Deviance = 594.3588 
#> GAMLSS-RS iteration 100: Global Deviance = 594.3538 
#> GAMLSS-RS iteration 101: Global Deviance = 594.3492 
#> GAMLSS-RS iteration 102: Global Deviance = 594.3448 
#> GAMLSS-RS iteration 103: Global Deviance = 594.3407 
#> GAMLSS-RS iteration 104: Global Deviance = 594.3367 
#> GAMLSS-RS iteration 105: Global Deviance = 594.333 
#> GAMLSS-RS iteration 106: Global Deviance = 594.3294 
#> GAMLSS-RS iteration 107: Global Deviance = 594.3261 
#> GAMLSS-RS iteration 108: Global Deviance = 594.323 
#> GAMLSS-RS iteration 109: Global Deviance = 594.3201 
#> GAMLSS-RS iteration 110: Global Deviance = 594.3174 
#> GAMLSS-RS iteration 111: Global Deviance = 594.3148 
#> GAMLSS-RS iteration 112: Global Deviance = 594.3124 
#> GAMLSS-RS iteration 113: Global Deviance = 594.3102 
#> GAMLSS-RS iteration 114: Global Deviance = 594.3078 
#> GAMLSS-RS iteration 115: Global Deviance = 594.3056 
#> GAMLSS-RS iteration 116: Global Deviance = 594.3035 
#> GAMLSS-RS iteration 117: Global Deviance = 594.3016 
#> GAMLSS-RS iteration 118: Global Deviance = 594.2998 
#> GAMLSS-RS iteration 119: Global Deviance = 594.2982 
#> GAMLSS-RS iteration 120: Global Deviance = 594.2966 
#> GAMLSS-RS iteration 121: Global Deviance = 594.2951 
#> GAMLSS-RS iteration 122: Global Deviance = 594.2938 
#> GAMLSS-RS iteration 123: Global Deviance = 594.2925 
#> GAMLSS-RS iteration 124: Global Deviance = 594.2913 
#> GAMLSS-RS iteration 125: Global Deviance = 594.2902 
#> GAMLSS-RS iteration 126: Global Deviance = 594.2892 
#> GAMLSS-RS iteration 127: Global Deviance = 594.2882 

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.95811    0.02332   41.09  < 2e-16 ***
#> x1           1.18156    0.04467   26.45  < 2e-16 ***
#> x2           0.22221    0.03637    6.11 2.01e-09 ***
#> ---
#> 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.9423     0.1264  15.370   <2e-16 ***
#> x3            1.3462     0.1583   8.503   <2e-16 ***
#> x4            1.4494     0.1667   8.697   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> No. of observations in the fit:  500 
#> Degrees of Freedom for the fit:  6
#>       Residual Deg. of Freedom:  494 
#>                       at cycle:  127 
#>  
#> Global Deviance:     594.2882 
#>             AIC:     606.2882 
#>             SBC:     631.5759 
#> ******************************************************************

# 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