bssn {bssn}R Documentation

Birnbaum-Saunders model based on Skew-Normal distribution


It provides the density, distribution function, quantile function, random number generator, likelihood function, moments and EM algorithm for Maximum Likelihood estimators for a given sample, all this for the three parameter Birnbaum-Saunders model based on Skew-Normal Distribution. Also, we have the random number generator for the mixture of Birbaum-Saunders model based on Skew-Normal distribution. Finally, the function mmmeth() is used to find the initial values for the parameters alpha and beta using modified-moment method.


dbssn(ti, alpha=0.5, beta=1, lambda=1.5)
pbssn(q,  alpha=0.5, beta=1, lambda=1.5)
qbssn(p,  alpha=0.5, beta=1, lambda=1.5)
rbssn(n,  alpha=0.5, beta=1, lambda=1.5)



vector of observations.


vector of quantiles.


vector of probabilities.


number of observations.


shape parameter.


scale parameter.


skewness parameter.


Are weights adding to 1. Each one of them (alpha, beta and lambda) must be a vector of length g if you want to generate a random numbers from a mixture distribution BSSN.


If alpha, sigma or lambda are not specified they assume the default values of 0.5, 1 and 1.5, respectively, belonging to the Birnbaum-Saunders model based on Skew-Normal distribution denoted by BSSN(0.5,1,1.5).

As discussed in Filidor et. al (2011) we say that a random variable T is distributed as an BSSN with shape parameter α>0, scale parameter β>0 and skewness parameter λ in R, if its probability density function (pdf) is given by

f(t)=2φ(a(t;α,β))Φ(λ a(t;α,β))A(t;α,β), t>0

where φ(.) and Φ(.) are the standard normal density and cumulative distribution function respectively. Also a(t;α,β)=(1/α)(√{t/β}-√{β/t}) and A(t;α,β)=t^{-3/2}(t+β)/(2α β^{1/2})


dbssn gives the density, pbssn gives the distribution function, qbssn gives the quantile function, rbssn generates a random sample and rmixbssn genrates a mixture random sample.

The length of the result is determined by n for rbssn, and is the maximum of the lengths of the numerical arguments for the other functions dbssn, pbssn and qbssn.


Rocio Maehara and Luis Benites


Vilca, Filidor; Santana, L. R.; Leiva, Victor; Balakrishnan, N. (2011). Estimation of extreme percentiles in Birnbaum Saunders distributions. Computational Statistics & Data Analysis (Print), 55, 1665-1678.

Santana, Lucia; Vilca, Filidor; Leiva, Victor (2011). Influence analysis in skew-Birnbaum Saunders regression models and applications. Journal of Applied Statistics, 38, 1633-1649.

See Also

EMbssn, momentsbssn, ozone, reliabilitybssn


## Not run: 
## Let's plot an Birnbaum-Saunders model based on Skew-Normal distribution!

## Density
 sseq <- seq(0,3,0.01)
 dens <- dbssn(sseq,alpha=0.2,beta=1,lambda=1.5)
 plot(sseq, dens,type="l", lwd=2,col="red", xlab="x", ylab="f(x)", main="BSSN Density function")

# Differing densities on a graph
# positive values of lambda
 y   <- seq(0,3,0.01)
 f1  <- dbssn(y,0.2,1,1)
 f2  <- dbssn(y,0.2,1,2)
 f3  <- dbssn(y,0.2,1,3)
 f4  <- dbssn(y,0.2,1,4)
 den <- cbind(f1,f2,f3,f4)

 matplot(y,den,type="l", col=c("deepskyblue4", "firebrick1", "darkmagenta", "aquamarine4"), ylab
 ="Density function",xlab="y",lwd=2,sub="(a)")

 legend(1.5,2.8,c("BSSN(0.2,1,1)", "BSSN(0.2,1,2)", "BSSN(0.2,1,3)","BSSN(0.2,1,4)"),
 col = c("deepskyblue4", "firebrick1", "darkmagenta", "aquamarine4"), lty=1:4,lwd=2,

#negative values of lambda
 y   <- seq(0,3,0.01)
 f1  <- dbssn(y,0.2,1,-1)
 f2  <- dbssn(y,0.2,1,-2)
 f3  <- dbssn(y,0.2,1,-3)
 f4  <- dbssn(y,0.2,1,-4)
 den <- cbind(f1,f2,f3,f4)

 matplot(y,den,type="l", col=c("deepskyblue4", "firebrick1", "darkmagenta", "aquamarine4"),
 ylab ="Density function",xlab="y",lwd=2,sub="(a)")
 legend(1.5,2.8,c("BSSN(0.2,1,-1)", "BSSN(0.2,1,-2)","BSSN(0.2,1,-3)", "BSSN(0.2,1,-4)"),
 col=c("deepskyblue4","firebrick1", "darkmagenta","aquamarine4"),lty=1:4,lwd=2,seg.len=2,

## Distribution Function
 sseq <- seq(0.1,6,0.05)
 df   <- pbssn(q=sseq,alpha=0.75,beta=1,lambda=3)
 plot(sseq, df, type = "l", lwd=2, col="blue", xlab="x", ylab="F(x)",
 main = "BSSN Distribution  function")

#Inverse Distribution Function
 prob <- seq(0,1,length.out = 1000)
 idf  <- qbssn(p=prob,alpha=0.75,beta=1,lambda=3)
 plot(prob, idf, type="l", lwd=2, col="gray30", xlab="x", ylab =
 expression(F^{-1}~(x)), mgp=c(2.3,1,.8))
 title(main="BSSN Inverse Distribution function")

#Random Sample Histogram
 sample <- rbssn(n=10000,alpha=0.75,beta=1,lambda=3)
 hist(sample,breaks = 70,freq = FALSE,main="")
 title(main="Histogram and True density")
 sseq   <- seq(0,8,0.01)
 dens   <- dbssn(sseq,alpha=0.75,beta=1,lambda=3)

##Random Sample Histogram for Mixture of BSSN
sample <- rmixbssn(n=1000,alpha,beta,lambda,pii)
hist(sample$y,breaks = 70,freq = FALSE,main="")
title(main="Histogram and True density")
temp   <- seq(min(sample$y), max(sample$y), length.out=1000)
lines(temp, (pii[1]*dbssn(temp, alpha[1], beta[1],lambda[1]))+(pii[2]*dbssn(temp, alpha[2]
, beta[2],lambda[2])), col="red", lty=3, lwd=3) # the theoretical density
lines(temp, pii[1]*dbssn(temp, alpha[1], beta[1],lambda[1]), col="blue", lty=2, lwd=3)
# the first component
lines(temp, pii[2]*dbssn(temp, alpha[2], beta[2],lambda[2]), col="green", lty=2, lwd=3)
# the second component

## End(Not run)

[Package bssn version 1.0 Index]