Nothing
pExtDep <- function(q, type, method="Parametric", model, par, plot=TRUE,
main, xlab, cex.lab, cex.axis, lwd,...){
# Checking quantile vector or matrix
if(is.vector(q)){
dim <- length(q)
q0 <- q
if(is.matrix(par)){
npar <- nrow(par)
p <- vector(length=npar)
if(any(q==0)){return(rep(NA,npar))}
}else if(is.vector(par)){
if(any(q==0)){return(NA)}
}
}else if(is.matrix(q)){
dim <- ncol(q)
qNA <- apply(q,1,function(x) any(x==0))
q0 <- q
nq0 <- nrow(q0)
q <- q[!qNA,]
if(is.matrix(q)){
nq <- nrow(q)
}else if(is.vector(q)){
nq <- 1
}
if(is.vector(par)){
p <- vector(length=nq)
}else if(is.matrix(par)){
p <- matrix(nrow=nrow(par), ncol=nq)
}
}else{
stop("q must be a vector or a matrix")
}
# Checking method
methods <- c("Parametric", "NonParametric")
if(!any(method == methods)){ stop("Wrong method specified")}
# Checking the type
types <- c("lower", "inv.lower", "upper")
if(method=="Parametric" && !any(type == types)){ stop("Wrong probability type specified")}
if(method=="Parametric"){
if(all(dim != c(2,3)) ){ stop("Dimension 2 or 3 only")}
# Checking the model
models <- c("HR", "ET", "EST")
if(!any(model == models)){ stop("model wrongly specified")}
# Checking model, par & dim are correctly specified
if(is.vector(par)){
if(!dim_ExtDep(model=model, par=par, dim=dim)){stop("Length of 'par' is incorrect")}
}else if(is.matrix(par)){
if(!dim_ExtDep(model=model, par=par[1,], dim=dim)){stop("Length of 'par' is incorrect")}
}else{
stop("'par' should be a vector or a matrix")
}
if(model == "HR"){
if(is.vector(q)){
if(is.vector(par)){
p <- p.hr(q=q, par=par, type=type)
}else if(is.matrix(par)){
p <- apply(par,1, p.hr, q=q, type=type)
}
}else if(is.matrix(q)){
if(is.vector(par)){
p <- apply(q,1, p.hr, par=par, type=type)
}else if(is.matrix(par)){
for(i in 1:nq){
p[,i] <- apply(par,1, p.hr, q=q[i,], type=type)
}
}
}
}else if(model == "ET"){
if(is.vector(q)){
if(is.vector(par)){
p <- p.et(q=q, par=par, type=type)
}else if(is.matrix(par)){
p <- apply(par,1, p.et, q=q, type=type)
}
}else if(is.matrix(q)){
if(is.vector(par)){
p <- apply(q,1, p.et, par=par, type=type)
}else if(is.matrix(par)){
for(i in 1:nq){
p[,i] <- apply(par,1, p.et, q=q[i,], type=type)
}
}
}
}else if(model == "EST"){
if(is.vector(q)){
if(is.vector(par)){
p <- p.est(q=q, par=par, type=type)
}else if(is.matrix(par)){
p <- apply(par,1, p.est, q=q, type=type)
}
}else if(is.matrix(q)){
if(is.vector(par)){
p <- apply(q,1, p.est, par=par, type=type)
}else if(is.matrix(par)){
for(i in 1:nq){
p[,i] <- apply(par,1, p.est, q=q[i,], type=type)
}
}
}
}
}else if(method == "NonParametric"){
if(is.vector(q)){
if(is.vector(par)){
p <- ph(w=q, beta=par)
}else if(is.matrix(par)){
for(i in 1:npar){
p[i] <- ph(w=q, beta=par[i,])
}
}
}else if(is.matrix(q)){
if(is.vector(par)){
for(i in 1:nq){
p[i] <- ph(w=q[i,], beta=par)
}
}else if(is.matrix(par)){
for(i in 1:nq){
for(j in 1:npar){
p[j,i] <- ph(w=q[i,], beta=par[j,])
}
}
}
}
}
if(is.matrix(par)){
if(plot){
if(missing(main) ){main <- ""}
if(missing(cex.lab)){cex.lab <- 1.4}
if(missing(cex.axis)){cex.axis <- 1.4}
if(missing(lwd)){lwd <- 2}
if(is.vector(p)){
if(missing(xlab)){
if(dim==2){
if(type=="lower"){
xlab <- paste("P(X<",q[1],",Y<",q[2], ")", sep="")
}else if(type=="inv.lower"){
xlab <- paste("1-P(X<",q[1],",Y<",q[2], ")", sep="")
}else if(type=="upper"){
xlab <- paste("P(X>",q[1],",Y>",q[2], ")", sep="")
}
}else if(dim==3){
if(type=="lower"){
xlab <- paste("P(X<",q[1],",Y<",q[2],",Z<",q[3], ")", sep="")
}else if(type=="inv.lower"){
xlab <- paste("1-P(X<",q[1],",Y<",q[2],",Z<",q[3], ")", sep="")
}else if(type=="upper"){
xlab <- paste("P(X>",q[1],",Y>",q[2],",Z>",q[3], ")", sep="")
}
}
}
Ke <- density(p) # KDE
Hi <- hist(p, prob=TRUE, col="lightgrey",
ylim=range(Ke$y), main=main, xlab=xlab ,
cex.lab=cex.lab, cex.axis=cex.axis, lwd=lwd, ...)
p_ic <- quantile(p, probs=c(0.025, 0.5, 0.975))
points(x=p_ic, y=c(0,0,0), pch=4, lwd=4)
points(x=mean(p), y=0, pch=16, lwd=4)
lines(Ke, lwd = 2, col = "dimgrey")
}else if(is.matrix(p)){
for(i in 1:nq){
if(missing(xlab)){
if(dim==2){
if(type=="lower"){
xlab <- paste("P(X<",q[i,1],",Y<",q[i,2], ")", sep="")
}else if(type=="inv.lower"){
xlab <- paste("1-P(X<",q[i,1],",Y<",q[i,2], ")", sep="")
}else if(type=="upper"){
xlab <- paste("P(X>",q[i,1],",Y>",q[i,2], ")", sep="")
}
}else if(dim==3){
if(type=="lower"){
xlab <- paste("P(X<",q[i,1],",Y<",q[i,2],",Z<",q[i,3], ")", sep="")
}else if(type=="inv.lower"){
xlab <- paste("1-P(X<",q[i,1],",Y<",q[i,2],",Z<",q[i,3], ")", sep="")
}else if(type=="upper"){
xlab <- paste("P(X>",q[i,1],",Y>",q[i,2],",Z>",q[i,3], ")", sep="")
}
}
}
Ke <- density(p[,i]) # KDE
Hi <- hist(p[,i], prob=TRUE, col="lightgrey",
ylim=range(Ke$y), main=main, xlab=xlab ,
cex.lab=cex.lab, cex.axis=cex.axis, lwd=lwd, ...)
p_ic <- quantile(p[,i], probs=c(0.025, 0.5, 0.975))
points(x=p_ic, y=c(0,0,0), pch=4, lwd=4)
points(x=mean(p[,i]), y=0, pch=16, lwd=4)
lines(Ke, lwd = 2, col = "dimgrey")
}
}
}
}
if(is.matrix(q0)){
if(sum(qNA) !=0){
p.temp <- p
if(is.vector(par)){
p <- vector(length=nq0)
p[!qNA] <- p.temp
p[qNA] <- rep(NA, sum(qNA))
}else if(is.matrix(par)){
p <- matrix(nrow=nrow(par), ncol=nq0)
p[,!qNA] <- p.temp
#p[,qNA] <- matrix(NA, nrow=nrow(par), ncol=sum(qNA))
}
}
}
return(p)
}
pFailure <- function(n, beta, u1, u2, mar1, mar2, type, plot, xlab, ylab, nlevels=10){
if(!(type %in% c("and", "or", "both")) ){stop(" 'type' must be 'and', 'or' or 'both'.")}
### Simulate data
w <- rh(n, beta)
r <- runif(n)^(-1)
y <- cbind(2*r*w, 2*r*(1-w))
### Compute ustar
u.star <- function(threshold, mar){
if(mar[3] == 0){
return( exp((threshold - mar[1])/mar[2]) )
}else{
return( (1 + mar[3]*(threshold-mar[1])/mar[2])^(1/mar[3]) )
}
}
ustar1 <- u.star(u1, mar1)
ustar2 <- u.star(u2, mar2)
ustar <- expand.grid(ustar1, ustar2)
if(type == "and"){
AND <- apply(ustar, 1, function(u) sum(y[,1] > u[1] & y[,2] > u[2]) )/n
AND <- matrix(AND, nrow=length(u1), ncol=length(u2), byrow=FALSE)
}
if(type == "or"){
OR <- apply(ustar, 1, function(u) sum(y[,1] > u[1] | y[,2] > u[2]) )/n
OR <- matrix(OR, nrow=length(u1), ncol=length(u2), byrow=FALSE)
}
if(type == "both"){
and_or <- function(u, y){
c(sum(y[,1] > u[1] & y[,2] > u[2]), sum(y[,1] > u[1] | y[,2] > u[2]))
}
AND_OR <- apply(ustar,1, and_or, y=y) / n
AND <- matrix(AND_OR[1,], nrow=length(u1), ncol=length(u2), byrow=FALSE)
OR <- matrix(AND_OR[2,], nrow=length(u1), ncol=length(u2), byrow=FALSE)
}
if(plot){
op <- par(mai = c(0.8, 0.8, 0.35, 0.1), mgp = c(2.5, 1, 0), cex.axis=2, cex.lab=2, cex.main=2)
on.exit(par(op))
if(type == "and"){
contour(u1, u2, AND, levels=round(seq(min(AND), max(AND), length=nlevels), 3),
xlab=xlab, ylab=ylab, main="P(X>x and Y>y)", labcex=1)
}
if(type == "or"){
contour(u1, u2, OR, levels=round(seq(min(OR), max(OR), length=nlevels), 3),
xlab=xlab, ylab=ylab, main="P(X>x or Y>y)", labcex=1)
}
if(type == "both"){
contour(u1, u2, AND, levels=round(seq(min(AND), max(AND), length=nlevels), 3),
xlab=xlab, ylab=ylab, main="P(X>x and Y>y)", labcex=1)
contour(u1, u2, OR, levels=round(seq(min(OR), max(OR), length=nlevels), 3),
xlab=xlab, ylab=ylab, main="P(X>x or Y>y)", labcex=1)
}
}
if(type == "and"){
return(list(AND=AND))
}
if(type == "or"){
return(list(OR=OR)) }
if(type == "both"){
return(list(AND=AND, OR=OR))
}
}
#########################################
#########################################
### Internal functions for pExtDep
#########################################
#########################################
## Husler-Reiss model
p.hr <- function(q, par, type){
p.hr.2d <- function(q, par, type){
I1 <- pnorm(par+log(q[2]/q[1])/(2*par))
I2 <- pnorm(par+log(q[1]/q[2])/(2*par))
if(type == "lower"){
return( exp(-1/q[1] * I1 - 1/q[2] * I2 ))
}else if(type == "inv.lower"){
return( 1 - exp(-1/q[1] * I1 - 1/q[2] * I2 ))
}else if(type == "upper"){
return( 1 - exp(-1/q[1]) - exp(-1/q[2]) + exp(-1/q[1] * I1 - 1/q[2] * I2 ) )
}
}
p.hr.3d <- function(q, par, type){
lambda12 <- par[1]; lambda13 <- par[2]; lambda23 <- par[3];
S1 <- matrix(c(4*lambda12^2, rep(2*(lambda12^2 + lambda13^2 - lambda23^2),2), 4*lambda13^2),nrow=2)
S2 <- matrix(c(4*lambda12^2, rep(2*(lambda12^2 + lambda23^2 - lambda13^2),2), 4*lambda23^2),nrow=2)
S3 <- matrix(c(4*lambda13^2, rep(2*(lambda13^2 + lambda23^2 - lambda12^2),2), 4*lambda23^2),nrow=2)
I1 <- pmesn(x=c(log(q[2]/q[1])+2*lambda12^2, log(q[3]/q[1])+2*lambda13^2),scale = S1 )
I2 <- pmesn(x=c(log(q[1]/q[2])+2*lambda12^2, log(q[3]/q[2])+2*lambda23^2),scale = S2 )
I3 <- pmesn(x=c(log(q[1]/q[3])+2*lambda13^2, log(q[2]/q[3])+2*lambda23^2),scale = S3 )
if(type == "lower"){
return( exp(-1/q[1] * I1 - 1/q[2] * I2 - 1/q[3] * I3 ))
}else if(type == "inv.lower"){
return( 1 - exp(-1/q[1] * I1 - 1/q[2] * I2 - 1/q[3] * I3 ))
}else if(type == "upper"){
return( 1 - exp(-1/q[1]) - exp(-1/q[2]) - exp(-1/q[3])
+ p.hr.2d(q[1:2], lambda12, type="lower") + p.hr.2d(q[c(1,3)], lambda13, type="lower")
+ p.hr.2d(q[2:3], lambda23, type="lower") - exp(-1/q[1] * I1 - 1/q[2] * I2 - 1/q[3] * I3 ))
}
}
if(any(par<=0)){stop('wrong value of parameters')}
dim <- length(q)
if(!dim_ExtDep(model="HR", par=par, dim=dim )){stop("Length of par is incorrect")}
if(dim == 2){
return( p.hr.2d(q=q, par=par, type=type))
}else if(dim == 3){
return( p.hr.3d(q=q, par=par, type=type))
}
}
## Extremal-t model
p.et <- function(q, par, type){
p.et.2d <- function(q, par, type){
if( ( (par[1] <= -1) || (par[1] >=1) || (par[2] <= 0))==1){stop("Extremal-t parameters ill defined")}
rho <- par[1]; nu <- par[2];
I1 <- pt(sqrt((nu+1)/(1-rho^2))*((q[2]/q[1])^(1/nu)-rho), df=nu+1)
I2 <- pt(sqrt((nu+1)/(1-rho^2))*((q[1]/q[2])^(1/nu)-rho), df=nu+1)
if(type == "lower"){
return( exp(-1/q[1] * I1 - 1/q[2] * I2 ))
}else if(type == "inv.lower"){
return( 1 - exp(-1/q[1] * I1 - 1/q[2] * I2 ))
}else if(type == "upper"){
return( 1/q[1] * (1-I1) + 1/q[2] * (1-I2) )
}
}
p.et.3d <- function(q, par, type){
if( (any(par[1:3] <= -1) || any(par[1:3] >=1) || (par[4] <= 0))==1){stop("Extremal-t parameters ill defined")}
rho12 <- par[1]; rho13 <- par[2]; rho23 <- par[3]; nu <- par[4];
R1 <- (rho23-rho12*rho13)/sqrt((1-rho12^2)*(1-rho13^2))
R2 <- (rho13-rho12*rho23)/sqrt((1-rho12^2)*(1-rho23^2))
R3 <- (rho12-rho13*rho23)/sqrt((1-rho13^2)*(1-rho23^2))
if(sum(eigen(matrix(c(1,R1,R1,1),ncol=2))$values<0)>=1){stop("negative definite partial correlation matrix")}
if(sum(eigen(matrix(c(1,R2,R2,1),ncol=2))$values<0)>=1){stop("negative definite partial correlation matrix")}
if(sum(eigen(matrix(c(1,R3,R3,1),ncol=2))$values<0)>=1){stop("negative definite partial correlation matrix")}
x1 <- sqrt((nu+1)/(1-rho12^2))*((q[2]/q[1])^(1/nu)-rho12)
y1 <- sqrt((nu+1)/(1-rho13^2))*((q[3]/q[1])^(1/nu)-rho13)
x2 <- sqrt((nu+1)/(1-rho12^2))*((q[1]/q[2])^(1/nu)-rho12)
y2 <- sqrt((nu+1)/(1-rho23^2))*((q[3]/q[2])^(1/nu)-rho23)
x3 <- sqrt((nu+1)/(1-rho13^2))*((q[1]/q[3])^(1/nu)-rho13)
y3 <- sqrt((nu+1)/(1-rho23^2))*((q[2]/q[3])^(1/nu)-rho23)
I1 <- pmest(x=c(x1,y1), scale=matrix(c(1,R1,R1,1),ncol=2), df=nu+1)
I2 <- pmest(x=c(x2,y2), scale=matrix(c(1,R2,R2,1),ncol=2), df=nu+1)
I3 <- pmest(x=c(x3,y3), scale=matrix(c(1,R3,R3,1),ncol=2), df=nu+1)
if(type == "lower"){
return( exp(-1/q[1] * I1 - 1/q[2] * I2 - 1/q[3] * I3 ))
}else if(type == "inv.lower"){
return( 1 - exp(-1/q[1] * I1 - 1/q[2] * I2 - 1/q[3] * I3 ))
}else if(type == "upper"){
return( 1 - exp(-1/q[1]) - exp(-1/q[2]) - exp(-1/q[3])
+ p.et.2d(q[1:2], c(rho12,nu), type="lower") + p.et.2d(q[c(1,3)], c(rho13,nu), type="lower")
+ p.et.2d(q[2:3], c(rho23,nu), type="lower") - exp(-1/q[1] * I1 - 1/q[2] * I2 - 1/q[3] * I3 ))
}
}
dim <- length(q)
if(!dim_ExtDep(model="ET", par=par, dim=dim )){stop("Length of par is incorrect")}
if(dim == 2){
return( p.et.2d(q=q, par=par, type=type))
}else if(dim == 3){
return( p.et.3d(q=q, par=par, type=type))
}
}
## Extremal Skew-t
pmextst <- function(x, scale=1, shape=rep(0,length(x)), df=1){
if(any(is.na(x)))
return(NA)
.C("pmextst", as.double(x), as.double(scale), as.double(df),
as.double(shape), out=double(1), NAOK=TRUE)$out
}
p.est <- function(q, par, type){
## Preliminary functions
Sigma <- function(rho){
p <- round(uniroot(function(x){length(rho)-choose(x,2)},lower=1, upper=10)$root)
Sig <- diag(p)
Sig[lower.tri(Sig)] = rho
Sig[row(Sig) < col(Sig)] = t(Sig)[row(Sig) < col(Sig)]
return(Sig)
}
Sigma_j <-function(rho,j){
Sig <- Sigma(rho)
return(Sig[-j,-j] - Sig[-j,j] %*% t(Sig[j,-j]) )
}
s_j <- function(rho,j){
p <- round(uniroot(function(x){length(rho)-choose(x,2)},lower=1, upper=10)$root)
k=p-1;
sigma_j <- Sigma_j(rho,j)
M <- diag(k)
diag(M) = sqrt(diag(sigma_j))
return( M )
}
Sigma_bar_j <- function(rho,j){
sigma_j <- Sigma_j(rho,j)
sj <- s_j(rho,j)
return( solve(sj) %*% t(sigma_j) %*% solve(sj) )
}
alpha_tilde <- function(alpha,j){
return(t(alpha[-j]))
}
alpha_bar_j <- function(rho,alpha,j){
Sig <- Sigma(rho)
sigma_j <- Sigma_j(rho,j)
Alpha_tilde <- alpha_tilde(alpha,j)
num <- alpha[j] + Sig[j,-j] %*% t(Alpha_tilde)
denom <- sqrt( 1 + Alpha_tilde %*% t(sigma_j) %*% t(Alpha_tilde) )
return(num/denom)
}
alpha_star_j <- function(rho,alpha,j){
Alpha_tilde <- alpha_tilde(alpha,j)
sj <- s_j(rho,j)
return( Alpha_tilde %*% sj )
}
tau_star_j <- function(rho,alpha,nu,j){
Sig <- Sigma(rho)
Alpha_tilde <- alpha_tilde(alpha,j)
return( sqrt(nu+1) * (alpha[j] + Sig[-j,j] %*% t(Alpha_tilde) ) )
}
nu_p <- function(rho,alpha,nu,j){
Alpha_bar <- alpha_bar_j(rho,alpha,j)
return( pt(sqrt(nu+1)*Alpha_bar,df=nu+1) * 2^(nu/2-1) * gamma((nu+1)/2) / sqrt(pi) )
}
x_bar <- function(x,rho,alpha,nu,j){
return( x * nu_p(rho,alpha,nu,j) )
}
p.est.2d <- function(q, par, type){
dim <- length(q)
if(dim != 2){stop("Dimension should be 2")}
rho <- par[1]
alpha <- par[2:3]
nu <- par[4]
if( ( (rho <= -1) || (rho >=1) || (nu <= 0))==1){stop("Extremal-t parameters ill defined")}
nu_1 <- nu_p(rho=rho, alpha=alpha, nu=nu, j=1)
nu_2 <- nu_p(rho=rho, alpha=alpha, nu=nu, j=2)
x1_bar <- q[1] * nu_1
x2_bar <- q[2] * nu_2
comp1_1 <- sqrt((nu+1)/(1-rho^2))*((x2_bar/x1_bar)^(1/nu)-rho)
comp1_2 <- sqrt((nu+1)/(1-rho^2))*((x1_bar/x2_bar)^(1/nu)-rho)
alpha_star_1 <- alpha_star_j(rho=rho, alpha=alpha, j=1)
alpha_star_2 <- alpha_star_j(rho=rho, alpha=alpha, j=2)
tau_star_1 <- tau_star_j(rho=rho, alpha=alpha, nu=nu, j=1)
tau_star_2 <- tau_star_j(rho=rho, alpha=alpha, nu=nu, j=2)
I1 <- pest(x=comp1_1,scale=1,shape=alpha_star_1,extended=tau_star_1,df=nu+1)
I2 <- pest(x=comp1_2,scale=1,shape=alpha_star_2,extended=tau_star_2,df=nu+1)
if(type == "lower"){
return( exp(-1/q[1] * I1 - 1/q[2] * I2 ))
}else if(type == "inv.lower"){
return( 1 - exp(-1/q[1] * I1 - 1/q[2] * I2 ))
}else if(type == "upper"){
return( 1/q[1] * (1-I1) + 1/q[2] * (1-I2) )
}
}
p.est.3d <- function(q, par, type){
dim <- length(q)
if(dim != 3){stop("Dimension should be 3")}
rho <- par[1:3]
alpha <- par[4:6]
nu <- par[7]
if( (any(rho <= -1) || any(rho >=1) || (nu <= 0))==1){stop("Extremal-t parameters ill defined")}
Sigma_bar_1 <- Sigma_bar_j(rho=rho, j=1)
Sigma_bar_2 <- Sigma_bar_j(rho=rho, j=2)
Sigma_bar_3 <- Sigma_bar_j(rho=rho, j=3)
if(sum(eigen(Sigma_bar_1)$values<0)>=1){stop("negative definite partial correlation matrix")}
if(sum(eigen(Sigma_bar_2)$values<0)>=1){stop("negative definite partial correlation matrix")}
if(sum(eigen(Sigma_bar_3)$values<0)>=1){stop("negative definite partial correlation matrix")}
nu_1 <- nu_p(rho=rho, alpha=alpha, nu=nu, j=1)
nu_2 <- nu_p(rho=rho, alpha=alpha, nu=nu, j=2)
nu_3 <- nu_p(rho=rho, alpha=alpha, nu=nu, j=3)
x1_bar <- q[1] * nu_1
x2_bar <- q[2] * nu_2
x3_bar <- q[3] * nu_3
comp1_1 <- sqrt((nu+1)/(1-rho[1]^2))*((x2_bar/x1_bar)^(1/nu)-rho[1])
comp2_1 <- sqrt((nu+1)/(1-rho[2]^2))*((x3_bar/x1_bar)^(1/nu)-rho[2])
comp1_2 <- sqrt((nu+1)/(1-rho[1]^2))*((x1_bar/x2_bar)^(1/nu)-rho[1])
comp2_2 <- sqrt((nu+1)/(1-rho[3]^2))*((x3_bar/x2_bar)^(1/nu)-rho[3])
comp1_3 <- sqrt((nu+1)/(1-rho[2]^2))*((x1_bar/x3_bar)^(1/nu)-rho[2])
comp2_3 <- sqrt((nu+1)/(1-rho[3]^2))*((x2_bar/x3_bar)^(1/nu)-rho[3])
alpha_star_1 <- alpha_star_j(rho=rho, alpha=alpha, j=1)
alpha_star_2 <- alpha_star_j(rho=rho, alpha=alpha, j=2)
alpha_star_3 <- alpha_star_j(rho=rho, alpha=alpha, j=3)
tau_star_1 <- tau_star_j(rho=rho, alpha=alpha, nu=nu, j=1)
tau_star_2 <- tau_star_j(rho=rho, alpha=alpha, nu=nu, j=2)
tau_star_3 <- tau_star_j(rho=rho, alpha=alpha, nu=nu, j=3)
I1 <- pmest(x=c(comp1_1,comp2_1),scale=Sigma_bar_1,shape=alpha_star_1,extended=tau_star_1,df=nu+1)
I2 <- pmest(x=c(comp1_2,comp2_2),scale=Sigma_bar_2,shape=alpha_star_2,extended=tau_star_2,df=nu+1)
I3 <- pmest(x=c(comp1_3,comp2_3),scale=Sigma_bar_3,shape=alpha_star_3,extended=tau_star_3,df=nu+1)
if(type == "lower"){
return( exp(-1/q[1] * I1 - 1/q[2] * I2 - 1/q[3] * I3 ))
}else if(type == "inv.lower"){
return( 1 - exp(-1/q[1] * I1 - 1/q[2] * I2 - 1/q[3] * I3 ))
}else if(type == "upper"){
return( 1 - exp(-1/q[1]) - exp(-1/q[2]) - exp(-1/q[3])
+ p.est.2d(q[1:2], c(rho[1], alpha[1:2],nu), type="lower") + p.est.2d(q[c(1,3)], c(rho[2], alpha[c(1,3)],nu), type="lower")
+ p.est.2d(q[2:3], c(rho[3], alpha[2:3],nu), type="lower") - exp(-1/q[1] * I1 - 1/q[2] * I2 - 1/q[3] * I3 ))
}
}
dim <- length(q)
if(!dim_ExtDep(model="EST", par=par, dim=dim )){stop("Length of par is incorrect")}
if(dim == 2){
return( p.est.2d(q=q, par=par, type=type))
}else if(dim == 3){
return( p.est.3d(q=q, par=par, type=type))
}
}
### Definition of the angular probability measure
ph <- function(w, beta)
{
k <- length(beta) - 1
j <- 1:k
# res <- diff(beta) * dbeta(w, j, k-j+1)
res <- 0.5 * (diff(beta) + 1/k) * dbeta(w, j, k-j+1)
return(sum(res))
}
#########################################
#########################################
### END Internal functions for pExtDep
#########################################
#########################################
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.