DTP R Package
The DTP R package implements the family of Double Two-Piece distributions. The family of univariate double two–piece distributions is obtained by using a density–based transformation of unimodal symmetric continuous distributions with a shape parameter. The resulting distributions contain five interpretable parameters that control the mode, as well as the scale and shape in each direction.
References.
Bayesian modelling of skewness and kurtosis with Two-Piece Scale and shape distributions
Data from: Exploring Relationships in Body Dimensions
rm(list=ls())
#library(devtools)
#install_github("FJRubio67/DTP")
library(DTP)
# Heinz et al. (2003) data. Biometric Measurements, variable "waist girth".
data <- c(66.5,61.5,61.2,78.0,70.5,66.5,58.0,74.5,73.5,90.5,66.5,61.0,69.5,62.0,63.4,71.0,59.0,74.0,60.5,75.0,65.0,82.2,70.2,85.0,
75.0,63.0,70.0,72.6,78.7,61.0,64.2,63.1,85.0,66.0,65.0,63.2,83.6,68.0,75.7,65.5,66.9,67.9,67.8,63.2,73.1,65.3,65.6,62.7,
81.8,72.1,68.6,80.0,64.9,61.7,64.0,58.6,68.3,69.5,70.5,69.8,69.0,81.0,71.0,73.6,96.3,65.0,62.8,83.0,67.6,65.4,61.6,63.2,
59.4,62.7,66.6,64.6,65.7,63.0,78.5,68.2,68.5,68.8,69.0,78.0,62.0,65.8,65.0,64.9,71.5,61.4,75.4,67.2,69.5,57.9,64.4,62.3,
68.7,68.0,66.0,71.2,64.8,78.0,60.7,75.0,68.6,67.1,68.0,71.0,76.0,70.3,69.5, 101.5,67.2,68.7,66.0,67.5,65.2,64.2,61.0,60.2,
62.6,72.7,58.5,67.5,66.1,68.0,64.0,90.1,64.1,72.4,71.0,69.1,66.1,59.5,71.5,65.0,71.5,63.0,70.0,70.4,64.3,66.0,70.5,62.0,
67.6,65.0,64.0,70.5,79.2,71.1,74.7,62.9,70.6,61.6,74.2,65.5,67.7,72.9,72.1,93.4,60.7,69.2,75.5,85.6,71.6,75.5,73.7,67.9,
71.2,69.3,75.6,63.8,83.6,84.9,66.0,68.3,63.8,73.6,76.3,58.8,78.0,64.5,65.8,62.5,77.5,70.0,59.4,68.3,77.9,68.3,72.0,62.7,
58.0,72.7,85.4,67.3,70.8,67.4,66.0,63.6,71.7,73.6,66.9,65.5,70.4,80.3,65.0,60.4,75.3,83.4,72.2,70.1,69.8,66.4,61.7,76.4,
63.8,70.2,88.2,75.8,81.5,67.4,78.3,65.9,62.4,65.6,96.2,62.5,80.5,86.1,71.3,79.4,72.3,67.4,94.2,67.0,69.5,66.4,67.6,66.6,
69.7,72.0,64.5,72.3,69.1,74.7,60.8,66.9,68.8,63.7,71.2,79.6,77.9,69.6,86.0,69.9,63.5,57.9,72.2,80.4)
summary(data)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 57.90 64.75 68.30 69.80 72.75 101.50n <- length(data)
#########################################################################################
# DTP-t fit
#########################################################################################
# - log-likelihood function
loglik1 <- function(par){
if(par[2]>0 & par[3]>-1 & par[3]<1 & par[4]>0 & par[5]>0) return(- sum(ddtp(data,par[1],par[2],par[3],par[4],par[5],param="eps",dt,log=T)))
else return(Inf)
}
# Optimisation
OPT1 <- optim(c(65,5.5,0,10,5),loglik1,control=list(maxit=10000),method="BFGS")
AIC1 <- 2*OPT1$value + 2*5
BIC1 <- 2*OPT1$value + log(n)*5
MLE1 <- OPT1$par
#########################################################################################
# TPSC-t fit
#########################################################################################
# - log-likelihood function
loglik2 <- function(par){
if(par[2]>0 & par[3]>-1 & par[3]<1 & par[4]>0) return(- sum(ddtp(data,par[1],par[2],par[3],par[4],par[4],param="eps",dt,log=T)))
else return(Inf)
}
# Optimisation
OPT2 <- optim(c(65,5.5,0,10),loglik2,control=list(maxit=10000),method="BFGS")
AIC2 <- 2*OPT2$value + 2*4
BIC2 <- 2*OPT2$value + log(n)*4
MLE2 <- OPT2$par
#########################################################################################
# TPSH-t fit
#########################################################################################
# - log-likelihood function
loglik3 <- function(par){
if(par[2]>0 & par[3]>0 & par[4]>0) return(- sum(ddtp(data,par[1],par[2],0,par[3],par[4],param="eps",dt,log=T)))
else return(Inf)
}
# Optimisation
OPT3 <- optim(c(65,5.5,10,5),loglik3,control=list(maxit=10000),method="BFGS")
AIC3 <- 2*OPT3$value + 2*5
BIC3 <- 2*OPT3$value + log(n)*5
MLE3 <- OPT3$par
#########################################################################################
# Model comparison
#########################################################################################
# AIC
c(AIC1,AIC2,AIC3)## [1] 1738.837 1738.885 1750.749# BIC
c(BIC1,BIC2,BIC3)## [1] 1756.640 1753.127 1768.553# Estimated densities
tempf1 <- function(x) ddtp(x,MLE1[1],MLE1[2],MLE1[3],MLE1[4],MLE1[5],param="eps",dt)
tempf2 <- function(x) ddtp(x,MLE2[1],MLE2[2],MLE2[3],MLE2[4],MLE2[4],param="eps",dt)
tempf3 <- function(x) ddtp(x,MLE3[1],MLE3[2],0,MLE3[3],MLE3[4],param="eps",dt)
hist(data, breaks = 20 ,probability=T,xlim=c(50,105),ylim=c(0,0.08))
curve(tempf1,50,105,add=T,lwd=2,lty=1)
curve(tempf2,50,105,add=T,lwd=2,lty=2)
curve(tempf3,50,105,add=T,lwd=2,lty=3)
box()
legend(90, 0.07, c("DTP","TPSC", "TPSH"),
text.col = "black", lty = c(1, 2, 3),
merge = TRUE, bg = "gray90")