The function UPHN() defines the Unit-Power Half-Normal
distribution, a two parameter
distribution, for a gamlss.family object to be used in GAMLSS fitting
using the function gamlss().
UPHN(mu.link = "log", sigma.link = "log")Returns a gamlss.family object which can be used
to fit a COMPO distribution
in the gamlss() function.
The UPHN distribution with parameters \(\mu\) and \(\sigma\) has a support in \((0, 1)\) and density given by
\(f(x| \mu, \sigma) = \frac{2\mu}{\sigma x^2} \phi(\frac{1-x}{\sigma x}) (2 \Phi(\frac{1-x}{\sigma x})-1)^{\mu-1}\)
for \(0 < x < 1\), \(\mu > 0\) and \(\sigma > 0\).
Santoro, K. I., Gómez, Y. M., Soto, D., & Barranco-Chamorro, I. (2024). Unit-Power Half-Normal Distribution Including Quantile Regression with Applications to Medical Data. Axioms, 13(9), 599.
# Example 1
# Generating random values with
# known mu and sigma
require(gamlss)
mu <- 1.5
sigma <- 4.0
y <- rUPHN(1000, mu, sigma)
mod1 <- gamlss(y~1, sigma.fo=~1, family=UPHN,
control=gamlss.control(n.cyc=5000, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = 1113.654
#> GAMLSS-RS iteration 2: Global Deviance = -1168.265
#> GAMLSS-RS iteration 3: Global Deviance = -1419.277
#> GAMLSS-RS iteration 4: Global Deviance = -1451.251
#> GAMLSS-RS iteration 5: Global Deviance = -1455.558
#> GAMLSS-RS iteration 6: Global Deviance = -1456.126
#> GAMLSS-RS iteration 7: Global Deviance = -1456.233
#> GAMLSS-RS iteration 8: Global Deviance = -1456.245
#> GAMLSS-RS iteration 9: Global Deviance = -1456.243
#> GAMLSS-RS iteration 10: Global Deviance = -1456.24
#> GAMLSS-RS iteration 11: Global Deviance = -1456.239
#> GAMLSS-RS iteration 12: Global Deviance = -1456.238
exp(coef(mod1, what="mu"))
#> (Intercept)
#> 1.465446
exp(coef(mod1, what="sigma"))
#> (Intercept)
#> 3.986704
# Example 2
# Generating random values under some model
# A function to simulate a data set with Y ~ UPHN
gendat <- function(n) {
x1 <- runif(n)
x2 <- runif(n)
mu <- exp(0.75 - 0.69 * x1) # Approx 1.5
sigma <- exp(0.5 - 0.64 * x2) # Approx 1.20
y <- rUPHN(n, mu, sigma)
data.frame(y=y, x1=x1, x2=x2)
}
dat <- gendat(n=2000)
mod2 <- gamlss(y~x1, sigma.fo=~x2, family=UPHN, data=dat,
control=gamlss.control(n.cyc=5000, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = -1753.08
#> GAMLSS-RS iteration 2: Global Deviance = -1791.485
#> GAMLSS-RS iteration 3: Global Deviance = -1797.5
#> GAMLSS-RS iteration 4: Global Deviance = -1798.647
#> GAMLSS-RS iteration 5: Global Deviance = -1798.852
#> GAMLSS-RS iteration 6: Global Deviance = -1798.889
#> GAMLSS-RS iteration 7: Global Deviance = -1798.897
#> GAMLSS-RS iteration 8: Global Deviance = -1798.9
#> GAMLSS-RS iteration 9: Global Deviance = -1798.9
summary(mod2)
#> Warning: summary: vcov has failed, option qr is used instead
#> ******************************************************************
#> Family: c("UPHN", "Unit-Power Half-Normal")
#>
#> Call: gamlss(formula = y ~ x1, sigma.formula = ~x2, family = UPHN,
#> data = dat, control = gamlss.control(n.cyc = 5000, trace = TRUE))
#>
#> Fitting method: RS()
#>
#> ------------------------------------------------------------------
#> Mu link function: log
#> Mu Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.78686 0.04589 17.14 <2e-16 ***
#> x1 -0.83090 0.07735 -10.74 <2e-16 ***
#> ---
#> 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) 0.51018 0.02658 19.20 <2e-16 ***
#> x2 -0.64457 0.04455 -14.47 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ------------------------------------------------------------------
#> No. of observations in the fit: 2000
#> Degrees of Freedom for the fit: 4
#> Residual Deg. of Freedom: 1996
#> at cycle: 9
#>
#> Global Deviance: -1798.9
#> AIC: -1790.9
#> SBC: -1768.497
#> ******************************************************************