ret<- list(dcoeff=kendr,dir=direction,Pxx=PX)
}
if (kendr< -mxc) {
mxc<- kendr
ret<- list(dcoeff=kendr,dir=direction,Pxx=PX)
}
if (outsd) {
vv<- Fvv+F.n(-v,-v,offset = 0, smoothing = "none")-
(Fuu+F.n(-u,-u,offset = 0, smoothing = "none"))*
(PY*(1.0-kendr)+kendr)-
PX*((Fyy+F.n(-y,-y,offset = 0, smoothing = "none"))*
(1.0-kendr)+2*kendr)
ret[["sd"]]<- sqrt(var(vv)/n)/(PX*(1.0-PY))
z<- qnorm((1+eps)/2)
ret[["conf"]]<- c(ret[[1]]-ret[[4]]*z,ret[[1]]+ret[[4]]*z)
}
}}
}
return(ret)
}
#' Multivariate Kendall's tau
#'
#' The function kendrm evaluates the multivariate Kendall regression coefficient.
#' It describes how well the response vector y can be fitted by a function of the regressor variables
#' which increases or decreases as the regressors increases.
#'
#' @usage kendtaum(x,y,direction=NULL,out=1,outsd=TRUE,eps=0.95)
#' @param x data matrix of regressor variables
#' @param y data matrix of the response vector
#' @param out value 1: full output, value 0: reduced output of one coefficient
#' that is largest in absolute value, value 2: only the largest coefficient
#' @param direction vector of length d (d is number of regressors),
#' value 1 refers to regressors leading to increasing y whenever this regressor increases,
#' value -1 refers to regressors leading to decreasing y whenever this regressor increases.
#' If direction=NULL, then all coefficients are computed.
#' @param outsd logical. If TRUE, the estimated standard deviation and confidence intervals
#' are evaluated.  If FALSE, the evaluation is suppressed
#' @param eps real value, confidence level.
#' @return A list of the coefficients for several directions with components:
#'   \describe{
#'     \item{dcoeff}{Kendall regression coefficient}
#'     \item{dir}{direction vector}
#'     \item{Pxx}{fraction of \eqn{X \le \breve{X}} where \eqn{\breve{X}} has the same distribution as \eqn{X}}
#'     \item{sd}{standard deviation}
#'     \item{conf}{confidence interval}
#'   }
#' @references Eckhard Liebscher (2026). Kendall regression coefficient revisited.
#' @examples
#' library(MASS)
#' data <- gilgais
#' kendr(data[,1:3],data[,4],out=1)
#' @export
kendtaum<-function(x,outsd=TRUE,eps=0.95){
n<- nrow(x)
d<- ncol(x)
x<- as.matrix(x)
if (n<4) {  stop("not enough sample items")}
if ((eps<=0.5)|(eps>=1)){ stop("wrong epsilon")}
if (d==1){stop("x has only one column")}
Fxx<- F.n(x,x,offset = 0, smoothing = "none")
d2<- 2^(d-1)
ret<- list(dcoeff=(2.0*d2*(sum(Fxx)-1.0)/(n-1.0))-1.0)/(d2-1.0)
if (outsd){
vv<- Fxx+F.n(-x,-x,offset = 0, smoothing = "none")
z<- qnorm((1+eps)/2)
ret[["sd"]]<- 2.0*d2*sqrt(var(vv)/n)/(d2-1.0)
ret[["conf"]]<- c(ret[[1]]-ret[[2]]*z,ret[[1]]+ret[[2]]*z)
}
return(ret)
}
library(MASS)
data<- gilgais
zetac(data[,1],data[,2],method="all")
psider<- function(x,id,par=-1){
h<- abs(x)
if (par<=0) {
if (id==3) {par<- 1.5}
if (id==4) {par<- 0.5}}
if (id==1){return(2.0*x)}   # Spearman-coefficient
if (id==2){return(ifelse(x<0,-1.0,ifelse(x>0,1.0,0.0)))}  # Spearman's footrule
if (id==3){return(ifelse(x<0,-par*h^(par-1),ifelse(x>0,par*h^(par-1),0.0)))}  # power coefficient
if (id==4){                # Huber function
h0<- ifelse(h<par,x,par*sign(x))
return(h0)}
}
library(MASS)
data<- gilgais
zetac(data[,1],data[,2],method="all")
psider<- function(x,id,par=-1){
h<- abs(x)
if (par<=0) {
if (id==3) {par<- 1.5}
if (id==4) {par<- 0.5}}
if (id==1){return(2.0*x)}   # Spearman-coefficient
if (id==2){return(ifelse(x<0,-1.0,ifelse(x>0,1.0,0.0)))}  # Spearman's footrule
if (id==3){return(ifelse(x<0,-par*h^(par-1),ifelse(x>0,par*h^(par-1),0.0)))}  # power coefficient
if (id==4){                # Huber function
h0<- ifelse(h<par,x,par*sign(x))
return(h0)}
}
zetac(data[,1],data[,2],method="all")
zetac<- function(x,y,method="Spearman",methodF=1,parH=0.5,parp=1.5,outsd=TRUE,eps=0.95){
n<- length(x)
if (n!=length(y)){ stop("number of sample items in x and y are different")} #check the parameters
if (n<4) {stop("not enough sample items")}
if ((!is.vector(x))|(!is.vector(y))){ stop("data format wrong")}
if ((eps<=0.5)|(eps>=1)){ stop("wrong epsilon")}
if (!methodF %in% 1:3) methodF<- 1
### compute normed difference of ranks
if (methodF==1) {
u<- rank(x)/n
v<- rank(y)/n
}else{
if (methodF==2){
u<- rank(x)/(n+1)
v<- rank(y)/(n+1)
} else {
u<- rank(x)/sqrt(n^2-1)
v<- rank(y)/sqrt(n^2-1)
}}
uu<- u-v
if (method=="all") {method<- c("Spearman","footrule","power","Huber")}
if (outsd){ z<- qnorm((1+eps)/2)}   #quantile
### Spearman coefficient
if ("Spearman" %in% method) {
h<- psi(uu,1)
s<- 1.0-sum(h)/(n*psif(1))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(uu,1)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sqrt(var(h+(h2/n))/n)/psif(1)     #standard deviation
outl<- list(Spearman=s,Spearman_se=sd,Spearman_conf=c(s-sd*z,s+sd*z))
}else{
outl<- list(Spearman=s)
}
} else {outl<- list()}
### Spearman's footrule
if ("footrule" %in% method) {
h<- psi(uu,2)
s<- 1.0-sum(h)/(n*psif(2))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(uu,2)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sd(h+(h2/n))/psif(2)     #standard deviation
outl<- append(outl,list(footrule=s,footrule_se=sd,footrule_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(footrule=s))
}}
### power function coefficient
if ("power" %in% method) {
h<- psi(uu,3,parp)
s<- 1.0-sum(h)/(n*psif(3,parp))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(uu,3,parp)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sd(h+(h2/n))/psif(3,parp)     #standard deviation
outl<- append(outl,list(power=s,power_se=sd,power_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(power=s))
}}
### Huber function coefficient
if ("Huber" %in% method) {
h<- psi(uu,4,parH)
s<- 1.0-sum(h)/(n*psif(4,parH))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(uu,4,parH)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sd(h+(h2/n))/psif(4,parH)     #standard deviation
outl<- append(outl,list(Huber=s,Huber_se=sd,Huber_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(Huber=s))
}}
return(outl)  ## output: list of computed coefficients,
}
zetac(data[,1],data[,2],method="all")
nrow(data)
zetac<- function(x,y,method="Spearman",methodF=1,parH=0.5,parp=1.5,outsd=TRUE,eps=0.95){
n<- length(x)
if (n!=length(y)){ stop("number of sample items in x and y are different")} #check the parameters
if (n<4) {stop("not enough sample items")}
if ((!is.vector(x))|(!is.vector(y))){ stop("data format wrong")}
if ((eps<=0.5)|(eps>=1)){ stop("wrong epsilon")}
if (!methodF %in% 1:3) methodF<- 1
### compute normed difference of ranks
if (methodF==1) {
u<- rank(x)/n
v<- rank(y)/n
}else{
if (methodF==2){
u<- rank(x)/(n+1)
v<- rank(y)/(n+1)
} else {
u<- rank(x)/sqrt(n^2-1)
v<- rank(y)/sqrt(n^2-1)
}}
uu<- u-v
if (method=="all") {method<- c("Spearman","footrule","power","Huber")}
if (outsd){ z<- qnorm((1+eps)/2)}   #quantile
### Spearman coefficient
if ("Spearman" %in% method) {
h<- psi(uu,1)
s<- 1.0-sum(h)/(n*psif(1))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(uu,1)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sqrt(var(h+(h2/n))/n)/psif(1)     #standard deviation
outl<- list(Spearman=s,Spearman_se=sd,Spearman_conf=c(s-sd*z,s+sd*z))
}else{
outl<- list(Spearman=s)
}
} else {outl<- list()}
### Spearman's footrule
if ("footrule" %in% method) {
h<- psi(uu,2)
s<- 1.0-sum(h)/(n*psif(2))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(uu,2)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sqrt(var(h+(h2/n))/n)/psif(2)     #standard deviation
outl<- append(outl,list(footrule=s,footrule_se=sd,footrule_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(footrule=s))
}}
### power function coefficient
if ("power" %in% method) {
h<- psi(uu,3,parp)
s<- 1.0-sum(h)/(n*psif(3,parp))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(uu,3,parp)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sqrt(var(h+(h2/n))/n)/psif(3,parp)     #standard deviation
outl<- append(outl,list(power=s,power_se=sd,power_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(power=s))
}}
### Huber function coefficient
if ("Huber" %in% method) {
h<- psi(uu,4,parH)
s<- 1.0-sum(h)/(n*psif(4,parH))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(uu,4,parH)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sqrt(var(h+(h2/n))/n)/psif(4,parH)     #standard deviation
outl<- append(outl,list(Huber=s,Huber_se=sd,Huber_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(Huber=s))
}}
return(outl)  ## output: list of computed coefficients,
}
zetac(data[,1],data[,2],method="all")
zetac(data[,1],data[,2],method="power",parp=2)
zetac(data[,1],data[,2],method="power",parp=5)
psider<- function(x,id,par=-1){
h<- abs(x)
if (par<=0) {
if (id==3) {par<- 1.5}
if (id==4) {par<- 0.5}}
if (id==1){return(2.0*x)}   # Spearman-coefficient
if (id==2){return(sign(x))}  # Spearman's footrule
if (id==3){return(sign(x)*par*h^(par-1))}  # power coefficient
if (id==4){                # Huber function
h0<- ifelse(h<par,x,par*sign(x))
return(h0)}
}
theta <- 2
d <- 2
n <- 1000
clay <- claytonCopula(theta, dim = d)
library(copula)
theta <- 2
d <- 2
n <- 1000
clay <- claytonCopula(theta, dim = d)
data <- rCopula(n, clay)
zetac(data[,1],data[,2],method="all",outsd=F)
options(digits=15)
zetac(data[,1],data[,2],method="all",outsd=F)
theta <- 2
d <- 2
n <- 1000
clay <- claytonCopula(theta, dim = d)
data <- rCopula(n, clay)
zetac(data[,1],data[,2],method="all",outsd=F)
theta <- 2
d <- 2
n <- 1000
clay <- claytonCopula(theta, dim = d)
data <- rCopula(n, clay)
zetac(data[,1],data[,2],method="all",parp=3,outsd=F)
nw<- 100000
n<- 500
vest<- array(dim=c(nw,2))
#tic("total")
for (jj in 1:nw){
xy<- rCopula(n, clay)
zeta<- zetac(xy[,1],xy[,2],method="Spearman",outsd=TRUE)
#print(kendr)
vest[jj,1]<- unlist(zeta$Spearman)
vest[jj,2]<- unlist(zeta$Spearman_se)
cat(jj,"\r")
}
mean(vest[1:nw,1])
sd(vest[1:nw,1])
mean(vest[1:nw,2])
sd(vest[1:nw,2])
zetac<- function(x,y,method="Spearman",methodF=1,parH=0.5,parp=1.5,outsd=TRUE,eps=0.95){
n<- length(x)
if (n!=length(y)){ stop("number of sample items in x and y are different")} #check the parameters
if (n<4) {stop("not enough sample items")}
if ((!is.vector(x))|(!is.vector(y))){ stop("data format wrong")}
if ((eps<=0.5)|(eps>=1)){ stop("wrong epsilon")}
if (!methodF %in% 1:3) methodF<- 1
### compute normed difference of ranks
if (methodF==1) {
u<- rank(x)/n
v<- rank(y)/n
}else{
if (methodF==2){
u<- rank(x)/(n+1)
v<- rank(y)/(n+1)
} else {
u<- rank(x)/sqrt(n^2-1)
v<- rank(y)/sqrt(n^2-1)
}}
uu<- u-v
if (method=="all") {method<- c("Spearman","footrule","power","Huber")}
if (outsd){ z<- qnorm((1+eps)/2)}   #quantile
### Spearman coefficient
if ("Spearman" %in% method) {
h<- psi(uu,1)
s<- 1.0-sum(h)/(n*psif(1))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(u[j]-v[j],1)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sqrt(var(h+(h2/n))/n)/psif(1)     #standard deviation
outl<- list(Spearman=s,Spearman_se=sd,Spearman_conf=c(s-sd*z,s+sd*z))
}else{
outl<- list(Spearman=s)
}
} else {outl<- list()}
### Spearman's footrule
if ("footrule" %in% method) {
h<- psi(uu,2)
s<- 1.0-sum(h)/(n*psif(2))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(u[j]-v[j],2)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sqrt(var(h+(h2/n))/n)/psif(2)     #standard deviation
outl<- append(outl,list(footrule=s,footrule_se=sd,footrule_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(footrule=s))
}}
### power function coefficient
if ("power" %in% method) {
h<- psi(uu,3,parp)
s<- 1.0-sum(h)/(n*psif(3,parp))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(u[j]-v[j],3,parp)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sqrt(var(h+(h2/n))/n)/psif(3,parp)     #standard deviation
outl<- append(outl,list(power=s,power_se=sd,power_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(power=s))
}}
### Huber function coefficient
if ("Huber" %in% method) {
h<- psi(uu,4,parH)
s<- 1.0-sum(h)/(n*psif(4,parH))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(u[j]-v[j],4,parH)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- 2.0*sqrt(var(h+(h2/n))/n)/psif(4,parH)     #standard deviation
outl<- append(outl,list(Huber=s,Huber_se=sd,Huber_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(Huber=s))
}}
return(outl)  ## output: list of computed coefficients,
}
nw<- 100000
n<- 500
vest<- array(dim=c(nw,2))
#tic("total")
for (jj in 1:nw){
xy<- rCopula(n, clay)
zeta<- zetac(xy[,1],xy[,2],method="Spearman",outsd=TRUE)
#print(kendr)
vest[jj,1]<- unlist(zeta$Spearman)
vest[jj,2]<- unlist(zeta$Spearman_se)
cat(jj,"\r")
}
mean(vest[1:nw,1])
sd(vest[1:nw,1])
mean(vest[1:nw,2])
sd(vest[1:nw,2])
zetac<- function(x,y,method="Spearman",methodF=1,parH=0.5,parp=1.5,outsd=TRUE,eps=0.95){
n<- length(x)
if (n!=length(y)){ stop("number of sample items in x and y are different")} #check the parameters
if (n<4) {stop("not enough sample items")}
if ((!is.vector(x))|(!is.vector(y))){ stop("data format wrong")}
if ((eps<=0.5)|(eps>=1)){ stop("wrong epsilon")}
if (!methodF %in% 1:3) methodF<- 1
### compute normed difference of ranks
if (methodF==1) {
u<- rank(x)/n
v<- rank(y)/n
}else{
if (methodF==2){
u<- rank(x)/(n+1)
v<- rank(y)/(n+1)
} else {
u<- rank(x)/sqrt(n^2-1)
v<- rank(y)/sqrt(n^2-1)
}}
uu<- u-v
if (method=="all") {method<- c("Spearman","footrule","power","Huber")}
if (outsd){ z<- qnorm((1+eps)/2)}   #quantile
### Spearman coefficient
if ("Spearman" %in% method) {
h<- psi(uu,1)
s<- 1.0-sum(h)/(n*psif(1))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(u[j]-v[j],1)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- sqrt(var(h+(h2/n))/n)/psif(1)     #standard deviation
outl<- list(Spearman=s,Spearman_se=sd,Spearman_conf=c(s-sd*z,s+sd*z))
}else{
outl<- list(Spearman=s)
}
} else {outl<- list()}
### Spearman's footrule
if ("footrule" %in% method) {
h<- psi(uu,2)
s<- 1.0-sum(h)/(n*psif(2))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(u[j]-v[j],2)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- sqrt(var(h+(h2/n))/n)/psif(2)     #standard deviation
outl<- append(outl,list(footrule=s,footrule_se=sd,footrule_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(footrule=s))
}}
### power function coefficient
if ("power" %in% method) {
h<- psi(uu,3,parp)
s<- 1.0-sum(h)/(n*psif(3,parp))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(u[j]-v[j],3,parp)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- sqrt(var(h+(h2/n))/n)/psif(3,parp)     #standard deviation
outl<- append(outl,list(power=s,power_se=sd,power_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(power=s))
}}
### Huber function coefficient
if ("Huber" %in% method) {
h<- psi(uu,4,parH)
s<- 1.0-sum(h)/(n*psif(4,parH))  #coefficient
if (outsd){
h2<- rep(0.0,n)
for (j in 1:n){
h2<- h2+psider(u[j]-v[j],4,parH)*((u<=u[j])-u[j]-(v<=v[j])+v[j])
}
sd<- sqrt(var(h+(h2/n))/n)/psif(4,parH)     #standard deviation
outl<- append(outl,list(Huber=s,Huber_se=sd,Huber_conf=c(s-sd*z,s+sd*z)))
}else{
outl<- append(outl,list(Huber=s))
}}
return(outl)  ## output: list of computed coefficients,
}
for (jj in 1:nw){
xy<- rCopula(n, clay)
zeta<- zetac(xy[,1],xy[,2],method="Spearman",outsd=TRUE)
#print(kendr)
vest[jj,1]<- unlist(zeta$Spearman)
vest[jj,2]<- unlist(zeta$Spearman_se)
cat(jj,"\r")
}
mean(vest[1:nw,1])
sd(vest[1:nw,1])
mean(vest[1:nw,2])
sd(vest[1:nw,2])
shiny::runApp('O:/Eckhard/stat/R/zuverl/zens_dstat')
