x<- cbind(x,-x)
for (k in ka:0){
if (ka>0){  #all directions
direction<- 2*as.numeric(intToBits(k)[1:d])-1   #vector of signs, first element always plus 1
if (k<ka)  { # direction does not consist solely of ones
sel<- (direction<0)
im<- (1:d) + d * sel   #select the columns of x
u<- x[,im]
}else{
u<- x[,1:d]
}
} else {  #  one direction
sel<- (direction<0)
im<- (1:d) + d * sel
u<- x[,im]
}
v<- cbind(u,y)
Fuu<- F.n(u,u,offset = 0, smoothing = "none")
Fvv<- F.n(v,v,offset = 0, smoothing = "none")
Fyy<- F.n(y,y,offset = 0, smoothing = "none")
PX<-(sum(Fuu)-1.0)/(n-1.0)
PXY<-(sum(Fvv)-1.0)/(n-1.0)
PY<-(sum(Fyy)-1.0)/(n-1.0)
if ((PX<=0.0)|(PY>=1.0)) {
if (out==1){
ret[[k+1]]<- list(dcoeff="not defined",dir=direction)
}else{  stop("Pxx=0, coefficient not defined")}
}else{
kendr<- (PXY-PX*PY)/(PX*(1.0-PY))    #Kendall coefficient
if (out==1){  # output of all coefficients
ret[[k+1]]<- list(dcoeff=kendr,dir=direction,Pxx=PX)
}else{  ## reduced output of the largest absolute value
akendr<- abs(kendr)
if (akendr>mxc) {
mxc<- akendr
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)
sd<- sqrt(var(vv)/n)/(PX*(1.0-PY))
z<- qnorm((1+eps)/2)
if (out==1){
ret[[k+1]][["sd"]]<-sd
ret[[k+1]][["conf"]]<- c(kendrc-sd*z,kendrc+sd*z)
}else{
ret[["sd"]]<- sd
ret[["conf"]]<- c(kendrc-sd*z,kendrc+sd*z)
}
} # end outsd
}
} # end k-loop
return(ret)
}
kendrm(x,y,direction=NULL,out=1,outsd=TRUE)
kendrm<-function(x,y,direction=NULL,out=0,outsd=TRUE,eps=0.95){
n<- nrow(x)
d<- ncol(x)
if (n!=nrow(y)){ stop("number of sample items in x and y are different")}
if (n<4) {  stop("not enough sample items")}
if ((eps<=0.5)|(eps>=1)){ stop("wrong epsilon")}
if (ncol(y)<2) { stop("too few columns in matrix y")}
x<- as.matrix(x)
y<- as.matrix(y)
mxc<- -10.0 #maximum coefficient
if (is.null(direction)) {
ka<-(2^d-1)
} else {
if (any(direction==0)|length(direction)!=d){stop("parameter direction not valid")}
ka<- 0
out<- 2}
ret<- list()
x<- cbind(x,-x)
for (k in ka:0){
if (ka>0){  #all directions
direction<- 2*as.numeric(intToBits(k)[1:d])-1   #vector of signs, first element always plus 1
if (k<ka)  { # direction does not consist solely of ones
sel<- (direction<0)
im<- (1:d) + d * sel   #select the columns of x
u<- x[,im]
}else{
u<- x[,1:d]
}
} else {  #  one direction
sel<- (direction<0)
im<- (1:d) + d * sel
u<- x[,im]
}
v<- cbind(u,y)
Fuu<- F.n(u,u,offset = 0, smoothing = "none")
Fvv<- F.n(v,v,offset = 0, smoothing = "none")
Fyy<- F.n(y,y,offset = 0, smoothing = "none")
PX<-(sum(Fuu)-1.0)/(n-1.0)
PXY<-(sum(Fvv)-1.0)/(n-1.0)
PY<-(sum(Fyy)-1.0)/(n-1.0)
if ((PX<=0.0)|(PY>=1.0)) {
if (out==1){
ret[[k+1]]<- list(dcoeff="not defined",dir=direction)
}else{  stop("Pxx=0, coefficient not defined")}
}else{
kendr<- (PXY-PX*PY)/(PX*(1.0-PY))    #Kendall coefficient
if (out==1){  # output of all coefficients
ret[[k+1]]<- list(dcoeff=kendr,dir=direction,Pxx=PX)
}else{  ## reduced output of the largest absolute value
akendr<- abs(kendr)
if (akendr>mxc) {
mxc<- akendr
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)
sd<- sqrt(var(vv)/n)/(PX*(1.0-PY))
z<- qnorm((1+eps)/2)
if (out==1){
ret[[k+1]][["sd"]]<-sd
ret[[k+1]][["conf"]]<- c(kendr-sd*z,kendr+sd*z)
}else{
ret[["sd"]]<- sd
ret[["conf"]]<- c(kendr-sd*z,kendr+sd*z)
}
} # end outsd
}
} # end k-loop
return(ret)
}
kendrm(x,y,direction=NULL,out=1,outsd=TRUE)
kendrm(x,y,direction=NULL,out=1,outsd=F)
kendrm(x,y,direction=NULL,out=0,outsd=T)
kendrm<-function(x,y,direction=NULL,out=0,outsd=TRUE,eps=0.95){
n<- nrow(x)
d<- ncol(x)
if (n!=nrow(y)){ stop("number of sample items in x and y are different")}
if (n<4) {  stop("not enough sample items")}
if ((eps<=0.5)|(eps>=1)){ stop("wrong epsilon")}
if (ncol(y)<2) { stop("too few columns in matrix y")}
x<- as.matrix(x)
y<- as.matrix(y)
mxc<- -10.0 #maximum coefficient
if (is.null(direction)) {
ka<-(2^d-1)
} else {
if (any(direction==0)|length(direction)!=d){stop("parameter direction not valid")}
ka<- 0
out<- 2}
ret<- list()
x<- cbind(x,-x)
for (k in ka:0){
if (ka>0){  #all directions
direction<- 2*as.numeric(intToBits(k)[1:d])-1   #vector of signs, first element always plus 1
if (k<ka)  { # direction does not consist solely of ones
sel<- (direction<0)
im<- (1:d) + d * sel   #select the columns of x
u<- x[,im]
}else{
u<- x[,1:d]
}
} else {  #  one direction
sel<- (direction<0)
im<- (1:d) + d * sel
u<- x[,im]
}
v<- cbind(u,y)
Fuu<- F.n(u,u,offset = 0, smoothing = "none")
Fvv<- F.n(v,v,offset = 0, smoothing = "none")
Fyy<- F.n(y,y,offset = 0, smoothing = "none")
PX<-(sum(Fuu)-1.0)/(n-1.0)
PXY<-(sum(Fvv)-1.0)/(n-1.0)
PY<-(sum(Fyy)-1.0)/(n-1.0)
if ((PX<=0.0)|(PY>=1.0)) {
if (out==1){
ret[[k+1]]<- list(dcoeff="not defined",dir=direction)
}else{  stop("Pxx=0, coefficient not defined")}
}else{
kendr<- (PXY-PX*PY)/(PX*(1.0-PY))    #Kendall coefficient
ff<- TRUE
if (out==1){  # output of all coefficients
ret[[k+1]]<- list(dcoeff=kendr,dir=direction,Pxx=PX)
}else{  ## reduced output of the largest absolute value
if (out==2) {
akendr<- kendr
} else { akendr<-abs(kendr)}
if (akendr>mxc) {
mxc<- akendr
ret<- list(dcoeff=kendr,dir=direction,Pxx=PX)
} else {  ff<- FALSE}
}
if ((outsd)&ff) {
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)
sd<- sqrt(var(vv)/n)/(PX*(1.0-PY))
z<- qnorm((1+eps)/2)
if (out==1){
ret[[k+1]][["sd"]]<-sd
ret[[k+1]][["conf"]]<- c(kendr-sd*z,kendr+sd*z)
}else{
ret[["sd"]]<- sd
ret[["conf"]]<- c(kendr-sd*z,kendr+sd*z)
}
} # end outsd
}
} # end k-loop
return(ret)
}
kendrm(x,y,direction=NULL,out=0,outsd=T)
kendrm(x,-y,direction=NULL,out=0,outsd=T)
kendrm(x,-y,direction=NULL,out=2,outsd=T)
kendrm(x,y,direction=NULL,out=2,outsd=T)
kendrm(x,-y,direction=c(-1,1,1,1),outsd=T)
theta <- 2
d <- 4
dimx<- 2
clay <- claytonCopula(theta, dim = d)
data <- rCopula(n, clay)
nw<- 100000
n<- 500
vest<- array(dim=c(nw,2))
for (jj in 1:nw){
xy<- rCopula(n, clay)
h<- kendrm(xy[,1:dimx],xy[,c(3,4)],direction=c(1,1),outsd=TRUE)
vest[jj,1]<- unlist(h$dcoeff)
vest[jj,2]<- unlist(h$sd)
cat(jj,"\r")
}
mean(vest[1:nw,1])
sd(vest[1:nw,1])
mean(vest[1:nw,2])
sd(vest[1:nw,2])
theta <- 2
d <- 3
dimx<- 2
clay <- claytonCopula(theta, dim = d)
data <- rCopula(n, clay)
nw<- 100000
n<- 500
vest<- array(dim=c(nw,2))
for (jj in 1:nw){
xx<- rCopula(n, clay)
h<- kendtaum(xx[,1:d],outsd=TRUE)
vest[jj,1]<- unlist(h$dcoeff)
vest[jj,2]<- unlist(h$sd)
cat(jj,"\r")
}
xx
Fxx<- F.n(x,x,offset = 0, smoothing = "none")
Fxx<- F.n(xx,xx,offset = 0, smoothing = "none")
d2<- 2^(d-1)
d
(2.0*d2*(sum(Fxx)-1.0)/(n-1.0))-1.0)/(d2-1.0)
ret<- list(dcoeff=((2.0*d2*(sum(Fxx)-1.0)/(n-1.0))-1.0)/(d2-1.0))
ret
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)
}
theta <- 2
d <- 3
dimx<- 2
clay <- claytonCopula(theta, dim = d)
data <- rCopula(n, clay)
nw<- 100000
n<- 500
vest<- array(dim=c(nw,2))
for (jj in 1:nw){
xx<- rCopula(n, clay)
h<- kendtaum(xx[,1:d],outsd=TRUE)
vest[jj,1]<- unlist(h$dcoeff)
vest[jj,2]<- unlist(h$sd)
cat(jj,"\r")
}
mean(vest[1:nw,1])
sd(vest[1:nw,1])
mean(vest[1:nw,2])
sd(vest[1:nw,2])
library(MASS)
data <- gilgais
kendtaum(data[,1:4],outsd=TRUE)
x1<- seq(0,1,by=0.01)
x<- cbind(x1,5*x1,3*x1)
kendtaum(x)
shiny::runApp('Eckhard/stat/R/zuverl/neu')
v <- c(10, 20)
c(a, b) <- as.list(v)
v <- c(10, 20)
list(a, b) <- as.list(v)
ew<-1
var<- 4
vv<- c(my=ew,v=var)
vv
my
shiny::runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
exists("get_nr_vert")
exists("get_nr_vert")
exists("Ht")
[1] FALSE
exists("get_nr_vert")
exists("get_moments")
source("global.R")
shiny::runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
par<- c(2,10)
tau0max<- tau0_max(nv,par)
tau0max<- tau0_max(3,par)
fakmaxtau
par[1]*par[2]
tau0_max(3,par)
hbeg<- step<- tau0max/300
tau_seq <- seq(hbeg, tau0max, length.out = 300)
tau_seq
kosten_seq <- sapply(tau_seq, function(t)
Kosten(t, nv, par, input$c_A, input$c_P))
Kosten <- function(tau0, nr_vert, par, c_A, c_P){
integral <- integrate(
function(x){ 1 - vertfkt(x, nr_vert, par)},
lower = 0,
upper = tau0
)$value
return(((c_A - c_P)*vertfkt(tau0, nr_vert, par)+c_P)/integral)
}
kosten_seq <- sapply(tau_seq, function(t)
Kosten(t, nv, par, input$c_A, input$c_P))
tau0
integrate(function(x) 1 - vertfkt(x, 2, c(2,10)), lower = 0, upper = 50)
integrate(function(x) 1 - vertfkt(x, 2, c(2,10)), lower = 0, upper = 50)$value
integrate(function(x) 1 - vertfkt(x, 3, c(2,10)), lower = 0, upper = 50)$value
integrate(function(x) 1 - vertfkt(x, "3", c(2,10)), lower = 0, upper = 50)$value
kosten_seq <- sapply(tau_seq, function(t)
+     Kosten(t, nv, par, input$c_A, input$c_P))
integral
kosten_seq <- sapply(tau_seq, function(t)
Kosten(t, 3, par, 1000, 100))
kosten_seq
for (j in 1:2){
tau_seq <- seq(hbeg, tau0max, length.out = 300)
# Kostenwerte berechnen
kosten_seq <- sapply(tau_seq, function(t)
Kosten(t, nv, par, input$c_A, input$c_P))
if (j==1){
kend<- kosten_seq[300]     #Endwert Kosten
grenze<- max(fakmaxkosten*min(kosten_seq),2*kend)
if (kosten_seq[1]<= grenze) { pos<- 0
}else{
pos <- max(which(kosten_seq > grenze))} #Position zum Abschneiden
hbeg<- tau_seq[pos+1]
}else{
pos<- which.min(kosten_seq) # naeherungsweises Minimum
if (pos==1){ interv<- c(tau_seq[1],tau_seq[2])}else{
if (pos==300){ interv<- c(tau_seq[299],tau_seq[300])}else{
interv<- c(tau_seq[pos-1],tau_seq[pos+1])}
}}
}
for (j in 1:2){
tau_seq <- seq(hbeg, tau0max, length.out = 300)
# Kostenwerte berechnen
kosten_seq <- sapply(tau_seq, function(t)
Kosten(t, nv, par, 1000,100))
if (j==1){
kend<- kosten_seq[300]     #Endwert Kosten
grenze<- max(fakmaxkosten*min(kosten_seq),2*kend)
if (kosten_seq[1]<= grenze) { pos<- 0
}else{
pos <- max(which(kosten_seq > grenze))} #Position zum Abschneiden
hbeg<- tau_seq[pos+1]
}else{
pos<- which.min(kosten_seq) # naeherungsweises Minimum
if (pos==1){ interv<- c(tau_seq[1],tau_seq[2])}else{
if (pos==300){ interv<- c(tau_seq[299],tau_seq[300])}else{
interv<- c(tau_seq[pos-1],tau_seq[pos+1])}
}}
}
j
tau0
kosten_seq
for (t in tau_seq)Kosten(t, nv, par,1000,100)
t
Kosten(tau_seq[1], nv, par,1000,100)
nv
nv<- 3
Kosten(tau_seq[1], nv, par,1000,100)
par
kosten_seq <- sapply(tau_seq, function(t)
Kosten(t, nv, par, input$c_A, input$c_P))
kosten_seq <- sapply(tau_seq, function(t)
Kosten(t, nv, par, 1000,100))
kosten_seq
kend<- kosten_seq[300]     #Endwert Kosten
grenze<- max(fakmaxkosten*min(kosten_seq),2*kend)
runApp('Eckhard/stat/R/zuverl/neu')
Kosten(6,3,c(2,10),1000,100)
nv<- 3
par<- c(2,10)
solve_tau0(nv, par, interv,1000,100)
nv<- 3
par<- c(2,10)
solve_tau0(nv, par, c(6,7),1000,100)
runApp('Eckhard/stat/R/zuverl/neu')
nv  <- 3
par <- c(2, 10)
c_A <- 1000
c_P <- 500
tau0max <- tau0_max(nv, par)
message("tau0max: ", tau0max)
tau_seq <- seq(tau0max/300, tau0max, length.out = 300)
kosten_seq <- sapply(tau_seq, function(t) Kosten(t, nv, par, c_A, c_P))
message("min kosten: ", min(kosten_seq, na.rm=TRUE))
message("pos minimum: ", which.min(kosten_seq))
message("NA in kosten: ", sum(is.na(kosten_seq)))
plot(tau_seq, kosten_seq, type="l")
tau_seq
kosten_seq
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
runApp('Eckhard/stat/R/zuverl/neu')
shiny::runApp('O:/Eckhard/stat/R/zuverl/prozessfaeh')
runApp('O:/Eckhard/stat/R/zuverl/prozessfaeh')
runApp('Eckhard/stat/R/zuverl/instandhaltg')
runApp('Eckhard/stat/R/zuverl/instandhaltg')
runApp('Eckhard/stat/R/zuverl/instandhaltg')
runApp('Eckhard/stat/R/zuverl/instandhaltg')
runApp('Eckhard/stat/R/zuverl/instandhaltg')
p<- 0.5
q<- 0.8
k <- floor(log(p) / log(q)) + 1
k
q^4
q^3
p<- 0.64
k <- floor(log(p) / log(q)) + 1
k
runApp('Eckhard/stat/R/zuverl/instandhaltg')
runApp('Eckhard/stat/R/zuverl/instandhaltg')
runApp('Eckhard/stat/R/zuverl/instandhaltg')
runApp('O:/Eckhard/stat/R/zuverl/deskr_stat')
runApp('O:/Eckhard/stat/R/zuverl/zens_dstat')
runApp('O:/Eckhard/stat/R/zuverl/zens_dstat')
runApp('O:/Eckhard/stat/R/zuverl/beschleun_lebensd')
runApp('O:/Eckhard/stat/R/zuverl/beschleun_lebensd')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
is.finite(NA)
is.finite(2.3)
is.finite(NaN)
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
runApp('O:/Eckhard/stat/R/zuverl/lebensdauervert')
