The Birnbaum-Saunders family - Santos-Neto et al. (2012) (P8 Based on the the first Tweedie)
Source:R/BS10.R
BS10.RdThe function BS10() defines the Birnbaum-Saunders distribution,
a two-parameter distribution, for a gamlss.family object
to be used in GAMLSS fitting using the function gamlss().
Value
Returns a gamlss.family object which can be used to fit a
BS10 distribution in the gamlss() function.
Details
The Birnbaum-Saunders distribution with parameters mu and sigma
(where mu represents \(\tau\) and sigma represents \(\omega\))
has density given by
\(f(x|\mu,\sigma) = \frac{1}{\sqrt{2\pi}} \exp\left( -\frac{\sigma}{\mu^2} \left[ \frac{2x}{\mu^2} + \frac{\mu^2}{2x} - 2 \right] \right) \frac{[2x + \mu^2]\sqrt{\sigma}}{2\mu^2 \sqrt{x^3}}\)
for \(x>0\), \(\mu>0\) and \(\sigma>0\). In this parameterization, \(E(X) = \frac{\mu^2}{2} + \frac{\mu^4}{8\sigma}\) and \(Var(X) = \frac{\mu^6}{8\sigma} + \frac{5\mu^8}{64\sigma^2}\).
References
Santos-Neto, M., Cysneiros, F. J. A., Leiva, V., & Ahmed, S. E. (2012). On new parameterizations of the Birnbaum-Saunders distribution. Pakistan Journal of Statistics, 28(1), 1-26.
Examples
# Example 1
# Generating some random values with
# known mu and sigma
set.seed(12345)
y <- rBS10(n=100, mu=0.75, sigma=5)
# Fitting the model
require(gamlss)
mod1 <- gamlss(y~1, sigma.fo=~1, family=BS10)
#> GAMLSS-RS iteration 1: Global Deviance = 63.8507
#> GAMLSS-RS iteration 2: Global Deviance = 32.2503
#> GAMLSS-RS iteration 3: Global Deviance = -19.5566
#> GAMLSS-RS iteration 4: Global Deviance = -108.0651
#> GAMLSS-RS iteration 5: Global Deviance = -216.8605
#> GAMLSS-RS iteration 6: Global Deviance = -243.7451
#> GAMLSS-RS iteration 7: Global Deviance = -244.1092
#> GAMLSS-RS iteration 8: Global Deviance = -244.1098
# Extracting the fitted values for mu and sigma
# using the inverse link function
exp(coef(mod1, what="mu"))
#> (Intercept)
#> 0.7506788
exp(coef(mod1, what="sigma"))
#> (Intercept)
#> 4.324429
# Example 2
# Generating random values for a regression model
# A function to simulate a data set with Y ~ BS10
gendat <- function(n) {
x1 <- runif(n)
x2 <- runif(n)
mu <- exp(0.5 - 1.6 * x1) # Aprox 0.75
sigma <- exp(2.2 - 1.2 * x2) # Aprox 5
y <- rBS10(n=n, mu=mu, sigma=sigma)
data.frame(y=y, x1=x1, x2=x2)
}
set.seed(123)
dat <- gendat(n=200)
mod2 <- gamlss(y~x1, sigma.fo=~x2,
family=BS10, data=dat)
#> GAMLSS-RS iteration 1: Global Deviance = 134.9485
#> GAMLSS-RS iteration 2: Global Deviance = -148.728
#> GAMLSS-RS iteration 3: Global Deviance = -359.1576
#> GAMLSS-RS iteration 4: Global Deviance = -525.3344
#> GAMLSS-RS iteration 5: Global Deviance = -542.83
#> GAMLSS-RS iteration 6: Global Deviance = -542.9624
#> GAMLSS-RS iteration 7: Global Deviance = -542.9627
summary(mod2)
#> Warning: summary: vcov has failed, option qr is used instead
#> ******************************************************************
#> Family: c("BS10", "Birnbaum-Saunders - Tenth parameterization")
#>
#> Call: gamlss(formula = y ~ x1, sigma.formula = ~x2, family = BS10,
#> data = dat)
#>
#> Fitting method: RS()
#>
#> ------------------------------------------------------------------
#> Mu link function: log
#> Mu Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.51286 0.02075 24.71 <2e-16 ***
#> x1 -1.61829 0.02730 -59.27 <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) 2.0661 0.1997 10.346 <2e-16 ***
#> x2 -0.9245 0.3564 -2.594 0.0102 *
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> ------------------------------------------------------------------
#> No. of observations in the fit: 200
#> Degrees of Freedom for the fit: 4
#> Residual Deg. of Freedom: 196
#> at cycle: 7
#>
#> Global Deviance: -542.9627
#> AIC: -534.9627
#> SBC: -521.7694
#> ******************************************************************