Nothing
# Density functions of SSMN Distributions
#require(moments)
#require(truncdist)
dssmn <- function(x, location=0, scale=1, shape=0, nu= 1, gama=1, dp=NULL, family="sn")
{
if(family=="sn")
{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
}
if(scale<=0){stop("Parameter scale must be positive")}
if(any(is.na(x))){x <- x[-which(is.na(x))]}
z <- (x-location)/sqrt(scale)
y <- 2*dnorm(z)*pnorm(shape*z)/sqrt(scale)
}
if(family=="stn")
{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
if (nu<=0){stop("Parameter nu must be positive")}
if(scale<=0){stop("Parameter scale must be positive")}
if (any(is.na(x))){ x <- x[-which(is.na(x))] }
z <- (x-location)/sqrt(scale)
cte <- gamma((nu+1)/2)/gamma(nu/2)/sqrt(nu*pi)
pdft <- cte*(1+z^2/nu)^(-(nu+1)/2)
y <- 2*pdft*pnorm(shape*z)/sqrt(scale)
}
if(family=="stn")
{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
if(nu<=0){stop("Parameter nu must be positive")}
if(scale<=0){stop("Parameter scale must be positive")}
if(any(is.na(x))){x <- x[-which(is.na(x))]}
z <- (x-location)/sqrt(scale)
n <- length(x)
cte <- nu/((2*pi*scale)^(1/2))
d <- z^2
fdps <- cte*gamma(nu+1/2)*pgamma(1,nu+1/2,scale=2/d)/((d/2)^(nu+1/2))
fdps[which(z==0)] <- cte/(nu+1/2)
y <- 2*fdps*pnorm(shape*z)
}
if(family=="stn")
{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
gama <- dp[5]
}
if(nu<=0|nu>=1){stop("Parameter nu must be between 0 and 1.0")}
if(gama<=0|gama>=1){stop("Parameter gama must be between 0 and 1.0")}
if(scale<=0){stop("Parameter scale must be positive")}
if(any(is.na(x))){x <- x[-which(is.na(x))]}
z <- (x-location)/sqrt(scale)
z2 <- z*sqrt(gama)
y <- 2*(nu*dnorm(x, mean=location, sd=sqrt(scale/gama))+(1-nu)*dnorm(x,mean=location,sd=sqrt(scale)))*pnorm(shape*z)
}
if(family=="sep")
{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
if(nu<=0.5|nu>1){stop("Parameter nu must be 0.5 < nu <= 1.0")}
if(scale<=0){stop("Parameter scale must be positive")}
if(any(is.na(x))){x<- x[-which(is.na(x))]}
nu1 <- 1/(2*nu)
z <- (x-location)/sqrt(scale)
d <- z^2
cte <- nu/(2^nu1)/gamma(nu1)/sqrt(scale)
pdfep <- cte*exp(-0.5*d^(nu))
y <- 2*pdfep*pnorm(shape*z)
}
return(list(y=y,family=family))
}
###############################################################################
# Cumulative Distributions functions of SSMN Distributions
pssmn <- function(q, location=0, scale=1, shape=0, nu= 1, gama=1, dp=NULL, family="sn")
{
x <- q
nrep <- length(x)
if(family=="sn")
{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
}else{
if(scale<=0){stop("Parameter scale must be positive")}
delta <- shape/sqrt(1+shape^2)
z <- (x-location)/sqrt(scale)
p <- numeric(length(x))
for (i in 1:length(x))
{
if(abs(x[i]) == Inf){ p[i] <- (1+sign(x[i]))/2
}else{
covi <- matrix(c(1,-delta,-delta,1),2,2)
Phi2 <- pmnorm(c(z[i],0),c(0,0),covi)
p[i] <- 2*Phi2
}
}
y <- pmax(0,pmin(1,p))
}
}
if(family=="stn")
{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}else{
if(nu<=0){stop("Parameter nu must be positive")}
if(scale<=0){stop("Parameter scale must be positive")}
delta <- shape/sqrt(1+shape^2)
U <- rgamma(nrep,shape=nu/2,scale=2/nu)
z <- (x-location)/sqrt(scale)
p <- numeric(length(x))
for(i in 1:length(x))
{
if(abs(x[i]) == Inf){ p[i] <- (1+sign(x[i]))/2
}else{
Phi2 <- numeric(nrep)
for (j in 1:nrep)
{
kappaj <- 1/U[j]
fj <- (1+shape^2*kappaj)/((1+shape^2)*kappaj)
covi <- kappaj*matrix(c(1,-delta,-delta,fj),2,2)
Phi2[j] <- pmnorm(c(z[i],0),c(0,0),covi)
}
p[i] <- 2*mean(Phi2)
}
}
y <- pmax(0,pmin(1,p))
}
}
if(family=="ssl")
{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}else{
if(nu<=0){stop("Parameter nu must be positive")}
if(scale<=0){stop("Parameter scale must be positive")}
delta <- shape/sqrt(1+shape^2)
V <- runif(nrep)
U <- V^(1/nu)
z <- (x-location)/sqrt(scale)
p <- numeric(length(x))
for (i in 1:length(x))
{
if(abs(x[i]) == Inf) {p[i] <- (1+sign(x[i]))/2
}else{
Phi2 <- numeric(nrep)
for (j in 1:nrep)
{
kappaj <- 1/U[j]
fj <- (1+shape^2*kappaj)/((1+shape^2)*kappaj)
covi <- kappaj*matrix(c(1,-delta,-delta,fj),2,2)
Phi2[j] <- pmnorm(c(z[i],0),c(0,0),covi)
}
p[i] <- 2*mean(Phi2)
}
}
y <- pmax(0,pmin(1,p))
}
}
if(family=="scn")
{
if(!is.null(dp))
{
if(!missing(shape)) {stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}else{
if(nu<=0){stop("Parameter nu must be positive")}
if (scale<=0){stop("Parameter scale must be positive")}
delta <- shape/sqrt(1+shape^2)
V <- rbinom(nrep,1,nu)
U <- gama*V+1-V
z <- (x-location)/sqrt(scale)
p <- numeric(length(x))
for (i in 1:length(x))
{
if(abs(x[i]) == Inf) {p[i] <- (1+sign(x[i]))/2
}else{
Phi2 <- numeric(nrep)
for (j in 1:nrep)
{
kappaj <- 1/U[j]
fj <- (1+shape^2*kappaj)/((1+shape^2)*kappaj)
covi <- kappaj*matrix(c(1,-delta,-delta,fj),2,2)
Phi2[j] <- pmnorm(c(z[i],0),c(0,0),covi)
}
p[i] <- 2*mean(Phi2)
}
}
y <- pmax(0,pmin(1,p))
}
}
if(family=="sep") #tem um problema no codigo p[i] <-
{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}else{
if(nu<=0.5|nu>1){stop("Parameter nu must be 0.5 < nu <= 1.0")}
if(scale<=0){stop("Parameter scale must be positive")
}
# z <- (x-location)/sqrt(scale)
p <- numeric(length(x))
for (i in 1:length(x)){
p[i] <-
if(abs(x[i]) == Inf) (1+sign(x[i]))/2
else{
integrate(dsep, -Inf, x[i], location=location, scale=scale, shape = shape, nu = nu)$value
}
}
pmax(0,pmin(1,p))
}
}
}
dsep <- function(x, location=0, scale=1, shape=0, nu=1, dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu<- dp[4]
}
if (nu<0.5|nu>1) {
stop("Parameter nu must be between 0.5 and 1.0")
}
if(scale<0) {
stop("Parameter scale must be positive")
}
if (any(is.na(x))) x<- x[-which(is.na(x))]
nu1<-1/(2*nu)
z<-(x-location)/sqrt(scale)
d<-z^2
cte<-nu/(2^nu1)/gamma(nu1)/sqrt(scale)
yt<-cte*exp(-0.5*d^(nu))
y<-2*yt*pnorm(shape*z)
return(y)
}
#pssmn(1, location=0, scale=1, shape=0, nu= 1, gama=1, dp=NULL, family="sn")
##############################################################################
# Quantile functions of SSMN Distributions #Revisar
qssmn <- function(p, location=0, scale=1, shape=0, nu= 1, gama=1, dp=NULL, family="sn")
{
if(family=="sn")
{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
}else{
if(scale<=0){stop("Parameter scale must be positive")}
na <- is.na(p) | (p < 0) | (p > 1)
zero <- (p == 0)
one <- (p == 1)
q <- numeric(length(p))
aux1 <- which(na | one | zero)
if(length(p)!=length(aux1))
{
aux2 <- c(1:length(p))
# aux2<- ifelse(length(aux1)!=0,aux2[-c(aux1)],aux2)
for (i in aux2)
{
q[i] <-
if (shape==0) qnorm(p[i])
else optimize(function(x) abs(psnn(x,0,1,shape)-p[i]),c((1.125 - 2.575*sign(shape))*sign(shape),(1.125 + 2.575*sign(shape))*sign(shape)))$minimum
}
}
}
q <- replace(q, na, NA)
q <- replace(q, zero, -Inf)
q <- replace(q, one, Inf)
return(location+sqrt(scale)*q)
}
qstn <- function (p,location=0, scale=1, shape=0, nu=30, dp=NULL) {
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
else
{
if (nu<=0) {
stop("Parameter nu must be positive")
}
if (scale<=0) {
stop("Parameter scale must be positive")
}
na <- is.na(p) | (p < 0) | (p > 1)
zero <- (p == 0)
one <- (p == 1)
q <- numeric(length(p))
aux1 <- which(na | one | zero)
if(length(p)!=length(aux1)) {
aux2 <- c(1:length(p))
#aux2<- ifelse(length(aux1)!=0,aux2[-c(aux1)],aux2)
for (i in aux2){
if (shape==0) q[i] <- qt(p[i],nu)
else { qti=optimize(function(x) abs(pstn(x,0,1,shape,nu)-p[i]),c((1.125 - 2.575*sign(shape))*sign(shape),(1.125 + 2.575*sign(shape))*sign(shape)))$minimum
q[i] <- qti}
}
}
}
q <- replace(q, na, NA)
q <- replace(q, zero, -Inf)
q <- replace(q, one, Inf)
return(location+sqrt(scale)*q)
}
qssl <- function (p,location=0, scale=1, shape=0, nu=30, dp=NULL) {
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
else
{
if (nu<=0) {
stop("Parameter nu must be positive")
}
if (scale<=0) {
stop("Parameter scale must be positive")
}
na <- is.na(p) | (p < 0) | (p > 1)
zero <- (p == 0)
one <- (p == 1)
q <- numeric(length(p))
aux1 <- which(na | one | zero)
if(length(p)!=length(aux1)) {
aux2 <- c(1:length(p))
# aux2<- ifelse(length(aux1)!=0,aux2[-c(aux1)],aux2)
for (i in aux2){
q[i] <-
if (shape==0) optimize(function(x) abs(pssl(x,0,1,shape,nu)-p[i]),c(-3,3))$minimum
else optimize(function(x) abs(pssl(x,0,1,shape,nu)-p[i]),c((1.125 - 2.575*sign(shape))*sign(shape),(1.125 + 2.575*sign(shape))*sign(shape)))$minimum
}
}
}
q <- replace(q, na, NA)
q <- replace(q, zero, -Inf)
q <- replace(q, one, Inf)
return(location+sqrt(scale)*q)
}
qscn <-function (p,location=0, scale=1, shape=0, nu=1, gama=1, dp=NULL) {
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
gama <- dp[5]
}
else
{
if (nu<=0|nu>=1) {
stop("Parameter nu must be between 0 and 1.0")
}
if(gama<=0|gama>=1) {
stop("Parameter gama must be between 0 and 1.0")
}
if (scale<=0) {
stop("Parameter scale must be positive")
}
na <- is.na(p) | (p < 0) | (p > 1)
zero <- (p == 0)
one <- (p == 1)
q <- numeric(length(p))
aux1 <- which(na | one | zero)
if(length(p)!=length(aux1)) {
aux2 <- c(1:length(p))
# aux2<- ifelse(length(aux1)!=0,aux2[-c(aux1)],aux2)
for (i in aux2){
q[i] <-
if (shape==0) optimize(function(x) abs(pscn(x,0,1,shape,nu,gama)-p[i]),c(-3,3))$minimum
else optimize(function(x) abs(pscn(x,0,1,shape,nu,gama)-p[i]),c((1.125 - 2.575*sign(shape))*sign(shape),(1.125 + 2.575*sign(shape))*sign(shape)))$minimum
}
}
}
q <- replace(q, na, NA)
q <- replace(q, zero, -Inf)
q <- replace(q, one, Inf)
return(location+sqrt(scale)*q)
}
qsep <-function (p,location=0, scale=1, shape=0, nu=1, dp=NULL) {
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
else
{
if (nu<=0.5|nu>1) {
stop("Parameter nu must be 0.5 < nu <= 1.0")
}
if (scale<=0) {
stop("Parameter scale must be positive")
}
na <- is.na(p) | (p < 0) | (p > 1)
zero <- (p == 0)
one <- (p == 1)
q <- numeric(length(p))
aux1 <- which(na | one | zero)
if(length(p)!=length(aux1)) {
aux2 <- c(1:length(p))
# aux2<- ifelse(length(aux1)!=0,aux2[-c(aux1)],aux2)
for (i in aux2){
q[i] <-
if (shape==0) optimize(function(x) abs(psep(x,0,1,shape,nu)-p[i]),c(-3,3))$minimum
else optimize(function(x) abs(psep(x,0,1,shape,nu)-p[i]),c((1.125 - 2.575*sign(shape))*sign(shape),(1.125 + 2.575*sign(shape))*sign(shape)))$minimum
}
}
}
q <- replace(q, na, NA)
q <- replace(q, zero, -Inf)
q <- replace(q, one, Inf)
return(location+sqrt(scale)*q)
}
}
## Auxiliares
psnn<- function (x, location=0, scale=1, shape=0, dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
}
else
{
if(scale<0) {
stop("Parameter scale must be positive")
}
# z <- (x-location)/sqrt(scale)
p <- numeric(length(x))
for (i in 1:length(x)){
p[i] <-
if(abs(x[i]) == Inf) (1+sign(x[i]))/2
else{
integrate(dsnn, -Inf, x[i], location=location, scale=scale, shape = shape)$value
}
}
pmax(0,pmin(1,p))
}
}
dsnn <- function (x, location=0, scale=1, shape=0, dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
}
if(scale<0) {
stop("Parameter scale must be positive")
}
if (any(is.na(x))) x<- x[-which(is.na(x))]
z<-(x-location)/sqrt(scale)
y<-2*dnorm(z)*pnorm(shape*z)/sqrt(scale)
return(y)
}
pstn<- function (x, location=0, scale=1, shape=0, nu=30, dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
else
{
if (nu<0) {
stop("Parameter nu must be positive")
}
if (scale<0) {
stop("Parameter scale must be positive")
}
# z <- (x-location)/sqrt(scale)
p <- numeric(length(x))
for (i in 1:length(x)){
p[i] <-
if(abs(x[i]) == Inf) (1+sign(x[i]))/2
else{
integrate(dstn, -Inf, x[i], location=location, scale=scale, shape = shape, nu = nu)$value
}
}
pmax(0,pmin(1,p))
}
}
dstn <- function(x, location=0, scale=1, shape=0, nu=30, dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
if (nu<0) {
stop("Parameter nu must be positive")
}
if(scale<0) {
stop("Parameter scale must be positive")
}
if (any(is.na(x))) x<- x[-which(is.na(x))]
z<-(x-location)/sqrt(scale)
cte<-gamma((nu+1)/2)/gamma(nu/2)/sqrt(nu*pi)
yt<-cte*(1+z^2/nu)^(-(nu+1)/2)
y<-2*yt*pnorm(shape*z)/sqrt(scale)
return(y)
}
pssl<- function(x, location=0, scale=1, shape=0, nu=30, dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu<- dp[4]
}
else
{
if (nu<0) {
stop("Parameter nu must be positive")
}
if(scale<0) {
stop("Parameter scale must be positive")
}
# z <- (x-location)/sqrt(scale)
p <- numeric(length(x))
for (i in 1:length(x)){
p[i] <-
if(abs(x[i]) == Inf) (1+sign(x[i]))/2
else{
integrate(dssl, -Inf, x[i], location=location, scale=scale, shape = shape, nu = nu)$value
}
}
pmax(0,pmin(1,p))
}
}
dssl <- function(x, location=0, scale=1, shape=0, nu=30, dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
if (nu<0) {
stop("Parameter nu must be positive")
}
if(scale<0) {
stop("Parameter scale must be positive")
}
if (any(is.na(x))) x<- x[-which(is.na(x))]
z<-(x-location)/sqrt(scale)
n<-length(x)
cte=nu/((2*pi*scale)^(1/2))
d=z^2
fdps=cte*gamma(nu+1/2)*pgamma(1,nu+1/2,scale=2/d)/((d/2)^(nu+1/2))
fsl=2*fdps*pnorm(shape*z)
return(fsl)
}
pscn <- function(x, location=0, scale=1, shape=0, nu=1, gama=1, dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape)) {
stop("You cannot set both component parameters and dp") }
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
gama <- dp[5]
}
else
{
if (nu<0|nu>1) {
stop("Parameter nu must be between 0.5 and 1.0")
}
if(gama<0|gama>1) {
stop("Parameter gama must be between 0.5 and 1.0")
}
if(scale<0) {
stop("Parameter scale must be positive")
}
# z <- (x-location)/sqrt(scale)
p <- numeric(length(x))
for (i in 1:length(x)){
p[i] <-
if(abs(x[i]) == Inf) (1+sign(x[i]))/2
else{
integrate(dscn, -Inf, x[i], location=location, scale=scale, shape = shape, nu = nu, gama=gama)$value
}
}
pmax(0,pmin(1,p))
}
}
dscn <- function(x, location=0, scale=1, shape=0, nu=1,gama=1,dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
gama <- dp[5]
}
if (nu<0|nu>1) {
stop("Parameter nu must be between 0.5 and 1.0")
}
if(gama<0|gama>1) {
stop("Parameter gama must be between 0.5 and 1.0")
}
if(scale<0) {
stop("Parameter scale must be positive")
}
if (any(is.na(x))) x<- x[-which(is.na(x))]
z<-(x-location)/sqrt(scale)
z2<-z*sqrt(gama)
y<-2*(nu*dnorm(x,mean=location,sd=sqrt(scale/gama))+(1-nu)*dnorm(x,mean=location,sd=sqrt(scale)))*pnorm(shape*z)
return(y)
}
psep <- function(x, location=0, scale=1, shape=0, nu=1, dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape)) {
stop("You cannot set both component parameters and dp") }
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
else
{
if (nu<0.5|nu>1) {
stop("Parameter nu must be between 0.5 and 1.0")
}
if(scale<0) {
stop("Parameter scale must be positive")
}
# z <- (x-location)/sqrt(scale)
p <- numeric(length(x))
for (i in 1:length(x)){
p[i] <-
if(abs(x[i]) == Inf) (1+sign(x[i]))/2
else{
integrate(dsep, -Inf, x[i], location=location, scale=scale, shape = shape, nu = nu)$value
}
}
pmax(0,pmin(1,p))
}
}
dsep <- function(x, location=0, scale=1, shape=0, nu=1, dp=NULL)
{
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu<- dp[4]
}
if (nu<0.5|nu>1) {
stop("Parameter nu must be between 0.5 and 1.0")
}
if(scale<0) {
stop("Parameter scale must be positive")
}
if (any(is.na(x))) x<- x[-which(is.na(x))]
nu1<-1/(2*nu)
z<-(x-location)/sqrt(scale)
d<-z^2
cte<-nu/(2^nu1)/gamma(nu1)/sqrt(scale)
yt<-cte*exp(-0.5*d^(nu))
y<-2*yt*pnorm(shape*z)
return(y)
}
################################################################################
# Random Generation of SSMN distributions
rssmn <- function(n, location=0, scale=1, shape=0, nu= 1, gama=1, dp=NULL, family="sn")
{
if(family == "sn")
{
if(!is.null(dp))
{
if(!missing(shape)) {stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
}
if(scale<=0) {stop("Parameter scale must be positive")}
delta <- shape/sqrt(1+shape^2)
u1 <- rnorm(n)
u2 <- rnorm(n)
y <- location+sqrt(scale)*(delta*abs(u1)+sqrt(1-delta^2)*u2)
}
if(family == "stn")
{
#rstn <- function(n, location, scale, shape, nu, dp=NULL)
#{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
if (nu<=0) {stop("Parameter nu must be positive")}
if(scale<=0){stop("Parameter scale must be positive")}
if (nu<1){warning('Nu < 1 can generate values tending to infinite',call. = FALSE)}
u <- rgamma(n,nu/2,nu/2)
ku <- 1/u
shape1 <- shape*sqrt(ku)
scale1 <- scale*ku
delta <- shape1/sqrt(1+shape1^2)
u1 <- rnorm(n)
u2 <- rnorm(n)
y <- location+sqrt(scale1)*(delta*abs(u1)+sqrt(1-delta^2)*u2)
#return(y)
#}
}
if(family == "ssl")
{
#rssl <- function(n, location, scale, shape, nu, dp)
#{
if(!is.null(dp))
{
if(!missing(shape)){stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
if(nu<=0){stop("Parameter nu must be positive")}
if(scale<0){stop("Parameter scale must be positive")}
v <- runif(n,0,1)
u <- v^(1/nu)
ku <- 1/u
shape1 <- shape*sqrt(ku)
scale1 <- scale*ku
delta <- shape1/sqrt(1+shape1^2)
u1 <- rnorm(n)
u2 <- rnorm(n)
y <- location+sqrt(scale1)*(delta*abs(u1)+sqrt(1-delta^2)*u2)
#return(y)
#}
}
if(family == "scn")
{
#rscn <- function(n, location, scale, shape, nu, gama, dp)
#{
if(!is.null(dp)) {
if(!missing(shape))
stop("You cannot set both component parameters and dp")
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
gama <- dp[5]
}
if (nu<0|nu>1) {
stop("Parameter nu must be between 0 and 1.0")
}
if(gama<=0|gama>1) {
stop("Parameter gama must be between 0 and 1.0")
}
if(scale<=0) {
stop("Parameter scale must be positive")
}
uu <- rbinom(n,1,nu)
u <- gama*uu+1-uu
ku <- 1/u
shape1 <- shape*sqrt(ku)
scale1 <- scale*ku
delta <- shape1/sqrt(1+shape1^2)
u1 <- rnorm(n)
u2 <- rnorm(n)
y <- location+sqrt(scale1)*(delta*abs(u1)+sqrt(1-delta^2)*u2)
# return(y)
#}
}
if(family == "sep")
{
#rsep <- function (n, location, scale, shape, nu, dp)
#{
if(!is.null(dp))
{
if(!missing(shape)) {stop("You cannot set both component parameters and dp")}
location <- dp[1]
scale <- dp[2]
shape <- dp[3]
nu <- dp[4]
}
if (nu<=0.5|nu>1){stop("Parameter nu must be 0.5 < nu <= 1.0")}
if(scale<=0) {stop("Parameter scale must be positive")}
t1 <- rgamma(n,shape=nu/2,scale=2)
R <- t1^(1/(2*nu))
u <- sign(rnorm(n))
yep <- u*R
w <- rnorm(n)
y <- location+sqrt(scale)*yep*sign(shape*yep-w)
#return(y)
#}
}
return(list(y=y,family=family))
}
#rssmn(n=20, location=0, scale=1, shape=0, nu= 1, gama=1, dp=NULL, family="sn")
###############################################################################
# Negative of Log-Likelihood of SSMN Regression Models
ssmn.logL <- function(theta,y,X, family="sn")
{
if(family=="sn") {
# Theta=[beta,sigma2,lambda]
if(missing(X)) {X<-matrix(1,n,1)}
p<-signif(ncol(X),digits=7)
beta<-signif(theta[1:p],digits=7)
sigma2<-signif(theta[1+p],digits=7)
sigma<-signif(sqrt(sigma2),digits=7)
lambda<-signif(theta[p+2],digits=7)
vero<-signif(2*dnorm(y,X%*%beta,sigma)*pnorm(lambda*(y-X%*%beta)/sigma),digits=7)
#logvero<-signif(log(vero),digits=7)
logvero=-sum(log(vero))
}
if(family=="stn") {
# Theta=[beta,sigma2,lambda,nu]
if(missing(X)) {X<-matrix(1,n,1)}
p<-ncol(X)
beta<-signif(theta[1:p],digits=7)
sigma2<-signif(theta[p+1],digits=7)
lambda<-signif(theta[p+2],digits=7)
nu<-signif(theta[p+3],digits=7)
z=(y-X%*%beta)/sqrt(sigma2)
d<-z^2
aux<-signif(lambda*z,digits=7)
aux1<-signif(pmax(aux,-37),digits=7)
cnugama<-signif(2*gamma((nu+1)/2)/(gamma(nu/2)*sqrt(nu*pi)),digits=7)
vero<-signif(cnugama/sqrt(sigma2)*(1+d/nu)^(-(nu+1)/2)*pnorm(aux1),digits=7)
logvero=-sum(log(vero))
}
if(family=="ssl") {
#% =Theta=[beta,sigma^2,lambda,nu]
n<-length(y)
if(missing(X)) {X<-matrix(1,n,1)}
p<-ncol(X)
beta<-theta[1:p]
sigma2<-theta[p+1]
lambda<-theta[p+2]
nu<-theta[p+3]
z=(y-X%*%beta)/sqrt(sigma2)
d<-z^2
aux<-signif(lambda*z,digits=7)
aux1<-signif(pmax(aux,-37),digits=7)
cte=nu/((2*pi*sigma2)^0.5)
fdps=cte*gamma(nu+0.5)*pgamma(1,nu+1/2,scale=2/d)/((d/2)^(nu+0.5))
fdps[which(z==0)]=cte/(nu+0.5)
vero=2*fdps*pnorm(aux1)
logvero=-sum(log(vero))
}
if(family=="scn") {
#% Theta=[beta,sigma2,lambda,nu,gama]
if(missing(X)) {X<-matrix(1,n,1)}
p<-ncol(X)
beta<-theta[1:p]
mu=X%*%beta
sigma2<-theta[p+1]
lambda<-theta[p+2]
nu<-theta[p+3]
gama<-theta[p+4]
z=(y-X%*%beta)/sqrt(sigma2)
aux<-signif(lambda*z,digits=7)
aux1<-signif(pmax(aux,-37),digits=7)
vero<-signif(2*(nu*dnorm(y,mu,sqrt(sigma2/gama))+(1-nu)*dnorm(y,mu,sqrt(sigma2)))*pnorm(aux1),digits=7)
logvero=-sum(log(vero))
}
if(family=="sep") {
#% Theta=[beta,sigma2,lambda,nu]
p<-ncol(X)
beta<-theta[1:p]
sigma2<-theta[p+1]
lambda<-theta[p+2]
nu<-theta[p+3]
z=(y-X%*%beta)/sqrt(sigma2)
d<-z^2
aux<-signif(lambda*z,digits=7)
aux1<-signif(pmax(aux,-37),digits=7)
cnu<-2*nu/(2^(1/(2*nu))*sqrt(sigma2)*gamma(1/(2*nu)))
vero<-cnu*exp(-0.5*d^nu)*pnorm(aux1)
logvero=-sum(log(vero))
}
return(logvero)
}
############################################################################
# Observed Fisher's Information Matrix
ssmn.info <- function(y,x,theta, family="sn")
{
if(family=="sn") {
# THETA: (beta,sigma2,lambda)
n=length(y)
if(missing(x)) {x<-matrix(1,n,1)}
p<-ncol(x)
beta<-theta[1:p]
sigma2<-as.numeric(theta[p+1])
lambda<-as.numeric(theta[p+2])
res<-y-x%*%beta
res<-as.vector(res)
sigma<-sqrt(sigma2)
aux<-lambda*res/sigma
aux1<-pmax(aux,-37)
Wphi<-dnorm(aux1)/pnorm(aux1)
d<-res^2/sigma2
######################### I2(beta,sigma2,lambda) ###########################
Qbeta<-t(res)%*%res
Wphi1<- -Wphi*(aux1+Wphi)
Qwbeta<- t(res)%*%diag(Wphi1)%*%res
I2<-matrix(0,p+2,p+2)
I2[1:p,1:p]<- -lambda^2/sigma2*t(x)%*%diag(Wphi1)%*%x
I2[1:p,p+1]<- -lambda/(2*sigma^3)*t(x)%*%Wphi-lambda^2/(2*sigma2^2)*t(x)%*%diag(Wphi1)%*%res
I2[1:p,p+2]<- 1/sigma*t(x)%*%Wphi+lambda/sigma2*t(x)%*%diag(Wphi1)%*%res
I2[p+1,p+1]<- -3*lambda/(4*sigma^5)*t(res)%*%Wphi-lambda^2/(4*sigma2^3)*Qwbeta
I2[p+1,p+2]<- 1/(2*sigma^3)*t(res)%*%Wphi+lambda/(2*sigma2^2)*Qwbeta
I2[p+2,p+2]<- -Qwbeta/sigma2
I2[p+2,p+1]<- I2[p+1,p+2]
I2[p+1,1:p]<- t(I2[1:p,p+1])
I2[p+2,1:p]<- t(I2[1:p,p+2])
######################## I1(beta,sigma2,lambda) ###########################
I1<-matrix(0,p+2,p+2)
I1[1:p,1:p]<-1/(sigma2)*t(x)%*%x
I1[1:p,p+1]<-1/(sigma2^2)*t(x)%*%res
I1[p+1,p+1]<- -n/(2*sigma2^2)+1/(sigma2^3)*t(res)%*%res
I1[p+1,1:p]<-t(I1[1:p,p+1])
########
Info <- I1+I2
}
if(family=="stn") {
# THETA: (beta,sigma2,lambda,nu)
n<-nrow(x)
if(missing(x)) {x<-matrix(1,n,1)}
p<-ncol(x)
beta<-theta[1:p]
sigma2<-as.numeric(theta[p+1])
lambda<-as.numeric(theta[p+2])
nu<-as.numeric(theta[p+3])
mu<-x%*%beta
res<-as.vector(y-mu)
sigma<-sqrt(sigma2)
aux<-lambda*res/sigma
aux1<-pmax(aux,-37)
Wphi<-dnorm(aux1)/pnorm(aux1)
d<-res^2/sigma2
######################### I2(beta,sigma2,lambda) ###########################
Qbeta<-t(res)%*%res
Wphi1<- as.vector(-Wphi*(aux1+Wphi))
Qwbeta<-t(res)%*%diag(Wphi1)%*%res
I2<-matrix(0,p+3,p+3)
I2[1:p,1:p]<- -lambda^2/sigma2*t(x)%*%diag(Wphi1)%*%x
I2[1:p,p+1]<- -lambda/(2*sigma^3)*t(x)%*%Wphi-lambda^2/(2*sigma2^2)*t(x)%*%diag(Wphi1)%*%res
I2[1:p,p+2]<- 1/sigma*t(x)%*%Wphi+lambda/sigma2*t(x)%*%diag(Wphi1)%*%res
I2[p+1,p+1]<- -3*lambda/(4*sigma^5)*t(res)%*%Wphi-lambda^2/(4*sigma2^3)*Qwbeta
I2[p+1,p+2]<- 1/(2*sigma^3)*t(res)%*%Wphi+lambda/(2*sigma2^2)*Qwbeta
I2[p+2,p+2]<- -Qwbeta/sigma2
I2[p+2,p+1]<- I2[p+1,p+2]
I2[p+1,1:p]<- t(I2[1:p,p+1])
I2[p+2,1:p]<- t(I2[1:p,p+2])
######################## I1(beta,sigma2,lambda,nu) ###########################
I1<-matrix(0,p+3,p+3)
V<-(1+d/nu)^(-1)
Qvbeta<- t(res)%*%diag(V)%*%res
hnu<-psigamma((nu+1)/2,1)-psigamma(nu/2,1)+2/(nu^2)
I1[1:p,1:p]<- -(nu+1)/(nu*sigma2)*t(x)%*%diag(V)%*%(2/nu*diag(V)%*%diag(d)-diag(n))%*%x
I1[1:p,p+1]<- -(nu+1)/(nu*sigma2^2)*t(x)%*%diag(res)%*%(1/nu*diag(V)%*%diag(d)-diag(n))%*%V
I1[1:p,p+3]<- -1/(nu^2*sigma2)*t(x)%*%diag(res)%*%((nu+1)/nu*diag(V)%*%diag(d)-diag(n))%*%V
I1[p+1,p+1]<- -n/(2*(sigma2^2))+(nu+1)/(nu*sigma2^3)*Qvbeta-(nu+1)/(2*(nu^2)*(sigma2^4))*t(res)%*%(diag(res^3))%*%diag(V)%*%V
I1[p+1,p+3]<- -0.5/((nu*sigma2)^2)*t(res)%*%diag(V)%*%((nu+1)/nu*diag(V)%*%diag(d)-diag(n))%*%res
I1[p+3,p+3]<- -n/4*hnu-0.5/(nu^2*sigma2)*Qvbeta-0.5/((nu^4)*sigma2)*t(res)%*%diag(V)%*%((nu+1)*diag(V)%*%diag(d)-nu*(nu+2)*diag(n))%*%res
I1[p+3,1:p]<-I1[1:p,p+3]
I1[p+1,1:p]<-I1[1:p,p+1]
I1[p+3,p+1]<-I1[p+1,p+3]
########
I<-I1+I2
return(I)
}
if(family=="ssl") {
# THETA: (beta,sigma2,lambda,nu)
n<-nrow(x)
if(missing(x)) {x<-matrix(1,n,1)}
p<-ncol(x)
beta<-theta[1:p]
sigma2<-as.numeric(theta[p+1])
lambda<-as.numeric(theta[p+2])
nu<-as.numeric(theta[p+3])
res<-as.vector(y-x%*%beta)
sigma<-sqrt(sigma2)
aux<-lambda*res/sigma
aux1<-pmax(aux,-37)
Wphi<-dnorm(aux1)/pnorm(aux1)
d<-res^2/sigma2
######################### I2(beta,sigma2,lambda) ###########################
Qbeta<-t(res)%*%res
Wphi1<- as.vector(-Wphi*(aux1+Wphi))
Qwbeta<-t(res)%*%diag(Wphi1)%*%res
I2<-matrix(0,p+3,p+3)
I2[1:p,1:p]<- -lambda^2/sigma2*t(x)%*%diag(Wphi1)%*%x
I2[1:p,p+1]<- -lambda/(2*sigma^3)*t(x)%*%Wphi-lambda^2/(2*sigma2^2)*t(x)%*%diag(Wphi1)%*%res
I2[1:p,p+2]<- 1/sigma*t(x)%*%Wphi+lambda/sigma2*t(x)%*%diag(Wphi1)%*%res
I2[p+1,p+1]<- -3*lambda/(4*sigma^5)*t(res)%*%Wphi-lambda^2/(4*sigma2^3)*Qwbeta
I2[p+1,p+2]<- 1/(2*sigma^3)*t(res)%*%Wphi+lambda/(2*sigma2^2)*Qwbeta
I2[p+2,p+2]<- -Qwbeta/sigma2
I2[p+2,p+1]<- I2[p+1,p+2]
I2[p+1,1:p]<- t(I2[1:p,p+1])
I2[p+2,1:p]<- t(I2[1:p,p+2])
######################## I1(beta,sigma2,lambda,nu) ###########################
I1<-matrix(0,p+3,p+3)
P11<-pgamma(1,nu+0.5,scale=2/d)
P13<-pgamma(1,nu+1.5,scale=2/d)
P15<-pgamma(1,nu+2.5,scale=2/d)
IG1<-gamma(nu+0.5)*P11/((d/2)^(nu+0.5))
IG3<-gamma(nu+1.5)*P13/((d/2)^(nu+1.5))
IG5<-gamma(nu+2.5)*P15/((d/2)^(nu+2.5))
rp<-5000
u1<-matrix(0,n,rp)
u2<-u1
for (j in 1:n) {
u1[j,]<-rtrunc(rp,"gamma",a=0,b=1,shape=nu+1/2,rate=d[j]/2)
u2[j,]<-rtrunc(rp,"gamma",a=0,b=1,shape=nu+3/2,rate=d[j]/2)
}
#E11=gamma(nu+0.5)*((d/2)^(-(nu+0.5)))*P11*rowMeans(log(u1))
#E13=gamma(nu+1.5)*((d/2)^(-(nu+1.5)))*P13*rowMeans(log(u2))
#E21=gamma(nu+2.5)*((d/2)^(-(nu+2.5)))*P11*rowMeans(log(u1)^2)
E11=P11*rowMeans(log(u1))
E13=P13*rowMeans(log(u2))
E21=P11*rowMeans(log(u1)^2)
I1[1:p,1:p]<- -1/sigma2*t(x)%*%diag(-IG3/IG1+d*(IG5-IG3^2/IG1)/IG1)%*%x
I1[1:p,p+1]<- -1/(2*sigma2^2)*t(x)%*%diag(res/IG1)%*%(-2*IG3+d*(IG5-IG3^2/IG1))
I1[1:p,p+3]<- -1/(sigma2)*t(x)%*%diag(res/(IG1^2))%*%(IG1*E13-IG3*E11)
I1[p+1,p+1]<- -1/(2*sigma2^2)*t(n+(d/IG1))%*%(-2*IG3+0.5*d*(IG5-IG3^2/IG1))
I1[p+1,p+3]<- -1/(2*sigma2)*t(d/(IG1^2))%*%(IG1*E13-IG3*E11)
I1[p+3,p+3]<- n/(nu^2)-t(IG1*E21-E11^2)%*%(IG1^(-2))
I1[p+1,1:p]<- t(I1[1:p,p+1])
I1[p+3,1:p]<- t(I1[1:p,p+3])
I1[p+3,p+1]<- I1[p+1,p+3]
########
I<-I1+I2
return(I)
}
if(family=="scn") {
# THETA: (beta,sigma2,lambda,nu,gama)
n<-nrow(x)
if(missing(x)) {x<-matrix(1,n,1)}
p<-ncol(x)
beta<-theta[1:p]
sigma2<-as.numeric(theta[p+1])
lambda<-as.numeric(theta[p+2])
nu<-as.numeric(theta[p+3])
gama<-as.numeric(theta[p+4])
res<-y-x%*%beta
sigma<-sqrt(sigma2)
aux<-lambda*res/sigma
aux1<-pmax(aux,-37)
Wphi<-dnorm(aux1)/pnorm(aux1)
d<-res^2/sigma2
######################### I2(beta,sigma2,lambda) ###########################
Qbeta<-t(res)%*%res
Wphi1<- as.vector(-Wphi*(aux1+Wphi))
Qwbeta<-t(res)%*%diag(Wphi1)%*%res
I2<-matrix(0,p+4,p+4)
I2[1:p,1:p]<- -lambda^2/sigma2*t(x)%*%diag(Wphi1)%*%x
I2[1:p,p+1]<- -lambda/(2*sigma^3)*t(x)%*%Wphi-lambda^2/(2*sigma2^2)*t(x)%*%diag(Wphi1)%*%res
I2[1:p,p+2]<- 1/sigma*t(x)%*%Wphi+lambda/sigma2*t(x)%*%diag(Wphi1)%*%res
I2[p+1,p+1]<- -3*lambda/(4*sigma^5)*t(res)%*%Wphi-lambda^2/(4*sigma2^3)*Qwbeta
I2[p+1,p+2]<- 1/(2*sigma^3)*t(res)%*%Wphi+lambda/(2*sigma2^2)*Qwbeta
I2[p+2,p+2]<- -Qwbeta/sigma2
I2[p+2,p+1]<- I2[p+1,p+2]
I2[p+1,1:p]<- t(I2[1:p,p+1])
I2[p+2,1:p]<- t(I2[1:p,p+2])
######################## I1(beta,sigma2,lambda,nu,gama) #######################
I1<-matrix(0,p+4,p+4)
sigma<-sqrt(sigma2)
fi<-dnorm(y,x%*%beta,sigma)
figama<-dnorm(y,x%*%beta,sqrt(sigma2/gama))
fs<-nu*figama+(1-nu)*fi
aux2<-as.vector(res*(nu*gama*figama+(1-nu)*fi))
dfsbeta<-1/sigma2*t(diag(aux2)%*%x)
dlbeta<-matrix(0,p,n)
for (i in 1:n) {
dlbeta[,i]<-dfsbeta[,i]/fs[i]
}
dfssigma<-1/(2*sigma2)*(nu*figama*(gama*d-1)+(1-nu)*fi*(d-1))
dfsnu<-figama-fi
dfsgama<-nu/(2*gama)*figama*(1-gama*d)
d2lbetabeta<-matrix(0,p,p)
for (i in 1:n) {
auxn<-nu*gama*figama[i]+(1-nu)*fi[i]
aux<-1/(sigma2*fs[i])*(d[i]*(nu*(gama^2)*figama[i]+(1-nu)*fi[i])-auxn-d[i]/fs[i]*auxn^2)*(as.vector(x[i,]))%*%t(as.vector(x[i,]))
d2lbetabeta<-d2lbetabeta+aux
}
d2lsigmabeta<-matrix(0,p,1)
for (i in 1:n) {
d2sigmabeta<-res[i]/(2*sigma2^2)*(nu*gama*figama[i]*(gama*d[i]-3)+(1-nu)*fi[i]*(d[i]-3))*as.vector(x[i,])
aux<-(d2sigmabeta*fs[i]-dfssigma[i]*dfsbeta[,i])/(fs[i]^2)
d2lsigmabeta<-d2lsigmabeta+aux
}
d2lnubeta<-matrix(0,p,1)
for (i in 1:n) {
d2nubeta<-1/sigma2*(gama*figama[i]-fi[i])*res[i]*as.vector(x[i,])
aux<--1/(fs[i]^2)*dfsnu[i]*dfsbeta[,i]+1/fs[i]*d2nubeta
d2lnubeta<-d2lnubeta+aux
}
d2lgamabeta<-matrix(0,p,1)
for (i in 1:n) {
d2gamabeta<-nu/(2*sigma2)*figama[i]*(3-gama*d[i])*res[i]*as.vector(x[i,])
aux<--1/(fs[i]^2)*dfsgama[i]*dfsbeta[,i]+1/fs[i]*d2gamabeta
d2lgamabeta<-d2lgamabeta+aux;
}
d2fssigmasigma=1/(2*sigma2^2)*(-2*nu*gama*figama*d-2*(1-nu)*fi*d+fs+nu/2*figama*(gama*d-1)^2+(1-nu)/2*fi*(d-1)^2);
Isigmasigma<-1/(fs^2)*(d2fssigmasigma*fs-dfssigma^2)
d2fsnusigma<-1/(2*sigma2)*(d*(gama*figama-fi)-figama+fi)
Inusigma<-1/(fs^2)*(d2fsnusigma*fs-dfssigma*dfsnu)
d2fsgamasigma<-nu/(4*sigma2)*figama*(4*d-gama*(d^2)-1/gama)
Igamasigma<-1/(fs^2)*(d2fsgamasigma*fs-dfssigma*dfsgama)
d2fsnunu<-0
Inunu<- -(dfsnu/fs)^2
d2fsnugama<-1/2*figama*(1/gama-d)
Inugama<-1/(fs^2)*(d2fsnugama*fs-dfsnu*dfsgama)
d2fsgamagama<-nu/4*figama*((1/gama-d)^2-2/(gama^2))
Igamagama<-1/(fs^2)*(d2fsgamagama*fs-dfsgama^2)
I1[1:p,1:p]<--d2lbetabeta
I1[1:p,p+1]<--d2lsigmabeta
I1[1:p,p+3]<--d2lnubeta
I1[1:p,p+4]<--d2lgamabeta
I1[p+1,p+1]<--sum(Isigmasigma)
I1[p+1,p+3]<--sum(Inusigma)
I1[p+1,p+4]<--sum(Igamasigma)
I1[p+3,p+3]<--sum(Inunu)
I1[p+3,p+4]<--sum(Inugama)
I1[p+4,p+4]<--sum(Igamagama)
I1[p+1,1:p]<- t(I1[1:p,p+1])
I1[p+3,1:p]<- t(I1[1:p,p+3])
I1[p+4,1:p]<- t(I1[1:p,p+4])
I1[p+3,p+1]<- I1[p+1,p+3]
I1[p+4,p+1]<- I1[p+1,p+4]
I1[p+4,p+3]<- I1[p+3,p+4]
########
I<-I1+I2
return(I)
}
if(family=="sep") {
# THETA: (beta,sigma2,lambda,nu)
n<-nrow(x)
if(missing(x)) {x<-matrix(1,n,1)}
p<-ncol(x)
beta<-theta[1:p]
sigma2<-as.numeric(theta[p+1])
lambda<-as.numeric(theta[p+2])
nu<-as.numeric(theta[p+3])
res<-as.vector(y-x%*%beta)
sigma<-sqrt(sigma2)
aux<-lambda*res/sigma
aux1<-pmax(aux,-37)
Wphi<-dnorm(aux1)/pnorm(aux1)
d<-res^2/sigma2
################### I2(beta,sigma2,lambda) ########################
Qbeta<-t(res)%*%res
Wphi1<- as.vector(-Wphi*(aux1+Wphi))
Qwbeta<-t(res)%*%diag(Wphi1)%*%res
I2<-matrix(0,p+3,p+3)
I2[1:p,1:p]<- -lambda^2/sigma2*t(x)%*%diag(Wphi1)%*%x
I2[1:p,p+1]<- -lambda/(2*sigma^3)*t(x)%*%Wphi-lambda^2/(2*sigma2^2)*t(x)%*%diag(Wphi1)%*%res
I2[1:p,p+2]<- 1/sigma*t(x)%*%Wphi+lambda/sigma2*t(x)%*%diag(Wphi1)%*%res
I2[p+1,p+1]<- -3*lambda/(4*sigma^5)*t(res)%*%Wphi-lambda^2/(4*sigma2^3)*Qwbeta
I2[p+1,p+2]<- 1/(2*sigma^3)*t(res)%*%Wphi+lambda/(2*sigma2^2)*Qwbeta
I2[p+2,p+2]<- -Qwbeta/sigma2
I2[p+2,p+1]<- I2[p+1,p+2]
I2[p+1,1:p]<- t(I2[1:p,p+1])
I2[p+2,1:p]<- t(I2[1:p,p+2])
######################## I1(beta,sigma2,lambda,nu) ###########################
I1<-matrix(0,p+3,p+3)
dnu<- d^(nu-1)
I1[1:p,1:p]<- nu*(2*nu-1)/sigma2*t(x)%*%diag(dnu)%*%x
I1[1:p,p+1]<- nu^2/(sigma2^2)*t(x)%*%diag(res)%*%(dnu)
I1[1:p,p+3]<- -1/(sigma2)*t(x)%*%diag(res*dnu)%*%(matrix(1,n,1)+nu*log(d))
I1[p+1,p+1]<- 1/(2*(sigma2^2))*(nu*(nu+1)*sum(d^nu)-n)
I1[p+1,p+3]<- -1/(2*sigma2)*t(d^nu)%*%(matrix(1,n,1)+nu*log(d))
nu1<- 1/(2*nu)
I1[p+3,p+3]<- n/(nu^2)*(1+(log(2)+psigamma(nu1))/nu+(nu1^2)*psigamma(nu1,1))+1/2*t(d^nu)%*%((log(d))^2)
I1[p+1,1:p]<- t(I1[1:p,p+1])
I1[p+3,1:p]<- t(I1[1:p,p+3])
I1[p+3,p+1]<- I1[p+1,p+3]
########
I<-I1+I2
return(I)
}
return(Info)
}
#ssmn.info(y,x,theta, family="sn")
gamtruncrnd<- function(a,b,t,n,m)
{
# Referencia: ANNE PHILIPPE. Simulation of right and left truncated gamma
# distributions by mixtures. Statistics and Computing (1997) 7, 173-181
# a: parameter of location
# t: right truncation, f_t varies from 0 to t
# [n,m]: sample sample
# If X ~ TG(a,b,t) then X/t ~ TG(a,bt,1)
# N: number of replications to ensure that the probability of acceptance algorithm ## acceptance-rejection P(N) > p.
tp=qnorm(0.99,0,1)
N=ceiling((tp+sqrt(tp^2+4*b))^2/4) # ceiling(X): integer greater than X
N=min(N,171)
K=matrix(1:N)
M=1/sum(b^(K-1)/gamma(K))
Y=matrix(0,n,m)
for (j in 1:n)
{
for (k in 1:m)
{
u=runif(1,0,1)
x=betamixt(a,b,N)
aux=((1-x)*b)^(K-1)/gamma(K)
rhox=1/(exp(b*x)*sum(aux))
Y[j,k]=x
while (u*M > rhox)
{
u=runif(1,0,1)
x=betamixt(a,b,N)
aux=((1-x)*b)^(K-1)/gamma(K)
rhox=1/(exp(b*x)*sum(aux))
Y[j,k]=x
}
}
}
return(Y)
}
betamixt<- function (a,b,N)
{
# gera de uma distribui?ao de misturas de betas
wb=matrix(1,N,1)
wt=wb
for (j in 2:N)
{
wb[j]=wb[j-1]*b/(a+j-1)
}
for (j in 2:N)
{
wt[j]=wt[j-1]+wb[j]
}
wt=rbind(0,wt/wt[N])
u=runif(1,0,1)
dif=wt-u
k=length(which(dif<0))
Y=rbeta(1,a,k)
return(Y)
}
################################################################################
# Negative of Log-likelihood of SMN regression models
ftn<-function(nu,y,X,theta) {
n<-length(y)
p=ncol(X)
beta=theta[1:p]
sigma2=theta[p+1]
erro=y-X%*%beta
aux=1+erro^2/(nu*sigma2)
logfy=n*log(gamma((nu+1)/2))-n/2*log(nu)-n*log(gamma(nu/2))-(nu+1)/2*sum(log(aux))
return(-logfy)
}
fsl<-function(nu,y,X,theta){
n<-length(y)
p=ncol(X)
beta=theta[1:p]
sigma2=theta[p+1]
cte=nu/sqrt(2*pi*sigma2)*gamma(nu+0.5)
d=(y-X%*%beta)^2/sigma2
fy=cte*pgamma(1,nu+0.5,scale=2/d)/((d/2)^(nu+0.5))
fy[which(d==0)]=cte/(nu+1/2)
return(-sum(log(fy)))
}
fcn <- function(nugama,y,X,theta)
{
n <- length(y)
p <- ncol(X)
beta <- theta[1:p]
mu <- X%*%beta
sigma2 <- theta[p+1]
nu <- nugama[1]
gama <- nugama[2]
fy <- nu*dnorm(y,mu,sqrt(sigma2/gama))+(1-nu)*dnorm(y,mu,sqrt(sigma2))
return(-sum(log(fy)))
}
fep <- function(nu,y,X,theta0)
{
n <- length(y)
p <- ncol(X)
beta <- theta0[1:p]
sigma2 <- theta0[p+1]
d <- (y-X%*%beta)^2/sigma2
cnu <- nu/(2^(1/(2*nu))*sqrt(sigma2)*gamma(1/(2*nu)))
fy <- cnu*exp(-0.5*d^nu)
return(-sum(log(fy)))
envelSCN <- function(Y,mu,Sigma,nu,gama){
n=length(Y)
p=1
d2=(Y-mu)^2/Sigma
d2s=sort(d2)
xq2 <- qchisq(ppoints(n), p)
Xsim<-matrix(0,200,n)
for(i in 1:200){
u1<-rchisq(n, p)/gama
u2<-rchisq(n, p)
u3<-runif(n)
id<-(u3<nu)
u2[id]<-u1[id]
Xsim[i,]=u2}
Xsim2<-apply(Xsim,1,sort)
d21<-rep(0,n)
d22<-rep(0,n)
for(i in 1:n){
d21[i] <- quantile(Xsim2[i,],0.025)
d22[i] <- quantile(Xsim2[i,],0.975)}
d2med <-apply(Xsim2,1,mean)
fy <- range(d2s,d21,d22)
plot(xq2,d2s,xlab = expression(paste("Theoretical ",chi[1]^2, " quantiles")),
ylab="Sample values and simulated envelope",pch=20,ylim=fy)
par(new=T)
plot(xq2,d21,type="l",ylim=fy,xlab="",ylab="")
par(new=T)
plot(xq2,d2med,type="l",ylim=fy,xlab="",ylab="",lty="dashed")
par(new=T)
plot(xq2,d22,type="l",ylim=fy,xlab="",ylab="")
}
}
################################################################################
# Negative of log-likelihood of SSMN regression models, estimating nu/gama
logL <- function(theta,y,X, family)
{
if(family=="sn")
{
#% Theta=[beta,sigma2,lambda]
p <- ncol(X)
beta <- theta[1:p]
sigma2 <- theta[p+1]
lambda <- theta[p+2]
z <- (y-X%*%beta)/sqrt(sigma2)
d <- z^2
aux <- signif(lambda*z,digits=7)
aux1 <- signif(pmax(aux,-37),digits=7)
cte <- 2/sqrt(2*pi*sigma2)
vero <- cte*exp(-0.5*d)*pnorm(aux1)
}
if(family=="stn")
{
# Theta=[beta,sigma2,lambda,nu]
if(missing(X)) {X<-matrix(1,n,1)}
p <- ncol(X)
beta <- signif(theta[1:p],digits=7)
sigma2 <- signif(theta[p+1],digits=7)
lambda <- signif(theta[p+2],digits=7)
nu <- optimize(ftn,c(1,30),y=y,X=X,theta=theta)$minimum
z <- (y-X%*%beta)/sqrt(sigma2)
d <- z^2
aux <- signif(lambda*z,digits=7)
aux1 <- signif(pmax(aux,-37),digits=7)
cnugama <- signif(2*gamma((nu+1)/2)/(gamma(nu/2)*sqrt(nu*pi)),digits=7)
vero <- signif(cnugama/sqrt(sigma2)*(1+d/nu)^(-(nu+1)/2)*pnorm(aux1),digits=7)
}
if(family=="ssl")
{
#% =Theta=[beta,sigma^2,lambda,nu]
n <- length(y)
if(missing(X)) {X<-matrix(1,n,1)}
p <- ncol(X)
beta <- theta[1:p]
sigma2 <- theta[p+1]
lambda <- theta[p+2]
nu <- optimize(fsl,c(0.1,20),y,X,theta)$minimum
z <- (y-X%*%beta)/sqrt(sigma2)
d <- z^2
aux <- signif(lambda*z,digits=7)
aux1 <- signif(pmax(aux,-37),digits=7)
cte <- nu/((2*pi*sigma2)^0.5)
fdps <- cte*gamma(nu+0.5)*pgamma(1,nu+0.5,scale=2/d)/((d/2)^(nu+0.5))
fdps[which(z==0)]=cte/(nu+0.5)
vero <- 2*fdps*pnorm(aux1)
}
if(family=="scn")
{
#% Theta=[beta,sigma2,lambda]
if(missing(X)) {X<-matrix(1,n,1)}
p <- ncol(X)
beta <- theta[1:p]
mu <- X%*%beta
sigma2 <- theta[p+1]
lambda <- theta[p+2]
L1 <- rbind(1e-6,1e-6)
L2 <- rbind(0.999999,0.999999)
nugama0 <- matrix(c(0.5,0.5),2,1)
nugama <- optim(nugama0,fcn,method='L-BFGS-B',lower=L1,upper=L2,control=list(maxit=400),y=y,X=X,theta=theta)$par
nu <- nugama[1]
gama <- nugama[2]
z <- (y-mu)/sqrt(sigma2)
aux <- signif(lambda*z,digits=7)
aux1 <- signif(pmax(aux,-37),digits=7)
vero <- signif(2*(nu*dnorm(y,mu,sqrt(sigma2/gama))+(1-nu)*dnorm(y,mu,sqrt(sigma2)))*pnorm(aux1),digits=7)
}
if(family=="sep")
{
#% Theta=[beta,sigma2,lambda]
p <- ncol(X)
beta <- theta[1:p]
sigma2 <- theta[p+1]
lambda <- theta[p+2]
nu <- optimize(fep,c(0.55,1),tol=1e-6,y,X,theta)$minimum
z <- (y-X%*%beta)/sqrt(sigma2)
d <- z^2
aux <- signif(lambda*z,digits=7)
aux1 <- signif(pmax(aux,-37),digits=7)
cnu <- 2*nu/(2^(1/(2*nu))*sqrt(sigma2)*gamma(1/(2*nu)))
vero <- cnu*exp(-0.5*d^nu)*pnorm(aux1)
}
return(-sum(log(vero)))
}
################################################################################
# Simulated envelope
envel <- function(y,X,theta,family="sn",alpha=0.05){
n=length(y)
p=ncol(X)
mu=X%*%theta[1:p]
d2=as.numeric((y-mu)^2/theta[p+1])
d2s=sort(d2)
d2s=t(d2s)
xq2 <- qchisq(ppoints(n), 1)
replic=200
Xsim<-matrix(0,replic,n)
if (family=="sn"){
for(i in 1:replic) Xsim[i,]<-rchisq(n, 1)
}
if (family=="stn") {
nu=theta[p+3]
for(i in 1:replic) Xsim[i,]<-rf(n, 1,nu)
}
if (family=="ssl"){
nu=theta[p+3]
for(i in 1:replic){
aux=GeraSlash(n,nu)
Xsim[i,]<-aux^2}
}
if (family=="scn"){
nu=theta[p+3]
gama=theta[p+4]
for(i in 1:replic){
u1<-rchisq(n, 1)/gama
u2<-rchisq(n, 1)
u3<-runif(n)
id<-(u3<nu)
u2[id]<-u1[id]
Xsim[i,]=u2
}
}
if (family=="sep"){
for(i in 1:replic){
aux=GeraEP(n,nu)
Xsim[i,]<-aux^2}
}
Xsim2<-apply(Xsim,1,sort)
d21<-rep(0,n)
d22<-rep(0,n)
for(i in 1:n){
d21[i] <- quantile(Xsim2[i,],alpha/2)
d22[i] <- quantile(Xsim2[i,],1-alpha/2)}
d2med <-apply(Xsim2,1,mean)
fy <- range(d2s,d21,d22)
plot(xq2,d2s,xlab = expression(paste("Theoretical ",chi[1]^2, " quantiles")),
ylab="Sample values and simulated envelope",pch=20,ylim=fy)
par(new=T)
plot(xq2,d21,type="l",ylim=fy,xlab="",ylab="")
par(new=T)
plot(xq2,d2med,type="l",ylim=fy,xlab="",ylab="",lty="dashed")
par(new=T)
plot(xq2,d22,type="l",ylim=fy,xlab="",ylab="")
}
#theta = c(fit.ssmn$beta,fit.ssmn$sigma2,fit.ssmn$lambda)
#envel(y,x1,theta,family="sn",alpha=0.05)
GeraSlash<- function(n,nu){
u1 <- runif(n)
u2 <- u1^(1/(nu))
u3<-rep(0,n)
for(j in 1:n){u3[j] <- rnorm(1, 0, sd=1/sqrt(u2[j]))}
return(u3)
}
GeraEP<- function(n, nu)
{
t=rgamma(n,0.5*nu,0.5)
R=t^(1/(2*nu))
u=sign(rnorm(n))
y=u*R
return(y)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.