The function PTRTE() defines the Poisson-transmuted record type exponential distribution,
a two-parameter discrete distribution, for a gamlss.family object to be used in GAMLSS fitting
using the function gamlss().
PTRTE(mu.link = "log", sigma.link = "logit")Returns a gamlss.family object that can be used to fit the Poisson-transmuted record type exponential
distribution using the gamlss() function.
The Poisson-transmuted record type exponential distribution with parameters \(\mu\) and \(\sigma\) has support \(x = 0,1,2,\dots\) and probability mass function given by
$$f(x | \mu, \sigma) = \frac{\mu(\sigma x\mu + 1 + \mu - \sigma)}{(1+\mu)^{x+2}}$$
Parameter restrictions: \(\mu > 0\) and \(0 < \sigma < 1\)).
Note: we renamed the original parameters \(\theta\) and \(p\) to \(\mu\) and \(\sigma\) respectively
to implement this distribution within the gamlss framework.
Erbayram, T., & Akdogan, Y. (2025). A new discrete model generated from mixed Poisson transmuted record type exponential distribution. Ricerca di Matematica, 74, 1225–1247.
# Example 1
# Generating some random values with known mu and sigma
# logit_inv function
logit_inv <- function(x) exp(x) / (exp(x)+1)
y <- rPTRTE(n=100, mu=0.2, sigma=0.5)
# Fitting the model
library(gamlss)
mod1 <- gamlss(y~1, family=PTRTE)
#> GAMLSS-RS iteration 1: Global Deviance = 611.5689
# Extracting the fitted values for mu and sigma
exp(coef(mod1, what="mu"))
#> (Intercept)
#> 0.2233241
logit_inv(coef(mod1, what="sigma"))
#> (Intercept)
#> 0.6749215
# Example 2
# Generating random values under some model
# A function to simulate a data set with Y ~ PTRTE
gendat <- function(n) {
x1 <- runif(n)
x2 <- runif(n)
mu <- exp(2 + 1 * x1) # 12 approximately
sigma <- logit_inv(2 - 2 * x2) # 0.73 approximately
y <- rPTRTE(n=n, mu=mu, sigma=sigma)
data.frame(y=y, x1=x1, x2=x2)
}
set.seed(1234)
dat <- gendat(n=1000)
# Fitting the model
mod2 <- NULL
mod2 <- gamlss(y~x1, sigma.fo=~x2, family=PTRTE, data=dat,
control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 900.7988
#> GAMLSS-RS iteration 2: Global Deviance = 897.0743
#> GAMLSS-RS iteration 3: Global Deviance = 897.0698
#> GAMLSS-RS iteration 4: Global Deviance = 897.0698
summary(mod2)
#> Warning: summary: vcov has failed, option qr is used instead
#> ******************************************************************
#> Family: c("PTRTE", "Poisson-Transmuted Record Type Exponential")
#>
#> Call: gamlss(formula = y ~ x1, sigma.formula = ~x2, family = PTRTE,
#> data = dat, 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) 1.7557 0.1705 10.297 < 2e-16 ***
#> x1 1.3392 0.3459 3.872 0.000115 ***
#> ---
#> 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) 126.0 618.0 0.204 0.839
#> x2 -196.7 959.0 -0.205 0.837
#>
#> ------------------------------------------------------------------
#> No. of observations in the fit: 1000
#> Degrees of Freedom for the fit: 4
#> Residual Deg. of Freedom: 996
#> at cycle: 4
#>
#> Global Deviance: 897.0698
#> AIC: 905.0698
#> SBC: 924.7008
#> ******************************************************************
# Example 3 (Second data set of the article)
# European corn-borer count data reported by McGuire et al. (1957).
# The observed and fitted frequencies are given in Table 11 of
# Erbayram and Akdogan (2025), where the P-TRTE distribution is
# illustrated using this data set.
values <- 0:5
freq <- c(188, 83, 36, 14, 2, 1)
y <- rep(x=values, times=freq)
mod3 <- gamlss(y~1, sigma.fo=~1, family=PTRTE(),
control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 711.8002
exp(coef(mod3, what="mu"))
#> (Intercept)
#> 2.874795
logit_inv(coef(mod3, what="sigma"))
#> (Intercept)
#> 0.8640032