}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")
expression(tau[0])
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')
nv  <- get_nr_vert("Weibull")   # dein Testfall
par <- get_par("Weibull", input) # ggf. manuell: c(2, 100)
t_seq <- seq(0, 200, length.out = 500)
S_t <- zuverlfkt(t_seq, nv, par, tau = 0)
for (t in t_seq)zuverlfkt(t, nv, par, tau = 10)
nr_vert
nv
get_nr_vert("Weibull")
setwd("o:/Eckhard/stat/R/zuverl/neu")
# globale Konstanten
maxc_A<- 1000   # maximale Ausfallkosten
fakmaxtau<- 10  # Faktor fuer maximales tau0
fakmaxkosten<- 20 # Faktor fuer maximale Kosten
# Packages
library(data.table)
library(ggplot2)
library(shinycssloaders)
library(stats)
library(shiny)
library(survival)
library(bslib)
# Kenngrößen ermitteln
get_moments <- function(nr_vert, par) {
switch(nr_vert,
"1" = {  # Weibull
beta <- par[1]
tau <- par[2]
g1 <- gamma(1 + 1/beta)
g2 <- gamma(1 + 2/beta)
g3 <- gamma(1 + 3/beta)
g4 <- gamma(1 + 4/beta)
ew   <- g1*tau
var  <- (g2 - g1^2)*(tau^2)
skew <- (g3 - 3*g2*g1 + 2*g1^3) / (g2 - g1^2)^(3/2)
kurt <- (g4 - 4*g3*g1 + 6*g2*g1^2 - 3*g1^4) / (g2 - g1^2)^2 - 3
},
"2" = {  # Exponential
ew  <- 1 / par[1]
var <- 1 / par[1]^2
skew <- 2
kurt <- 6
},
"3" = {  # Gamma
ew  <- par[1] * par[2]
var <- par[1] * par[2]^2
skew <- 2 / sqrt(par[1])
kurt <- 6/par[1]
},
"4" = {  # LogNormal
ew  <- exp(par[1] + par[2]^2 / 2)
var <- (exp(par[2]^2) - 1) * exp(2 * par[1] + par[2]^2)
skew <- (exp(par[2]^2) + 2) * sqrt(exp(par[2]^2) - 1)
kurt <- exp(4*par[2]^2) + 2*exp(3*par[2]^2) + 3*exp(2*par[2]^2) - 6
},
"5" = {  # Weibull 3 Par.
beta <- par[1]
tau <- par[2]
g1 <- gamma(1 + 1/beta)
g2 <- gamma(1 + 2/beta)
g3 <- gamma(1 + 3/beta)
g4 <- gamma(1 + 4/beta)
ew   <- g1*tau+par[3]
var  <- (g2 - g1^2)*(tau^2)
skew <- (g3 - 3*g2*g1 + 2*g1^3) / (g2 - g1^2)^(3/2)
kurt <- (g4 - 4*g3*g1 + 6*g2*g1^2 - 3*g1^4) / (g2 - g1^2)^2 - 3
}
)
return(c(ew, var, skew, kurt))
}
#Nummer Verteilung
get_nr_vert<- function(vchar){
return(
switch(vchar,
"Weibullverteilung"      = 1,
"Exponentialverteilung"  = 2,
"Gammaverteilung"        = 3,
"LogNormalverteilung"    = 4,
"Weibullverteilung 3 Par." = 5
))
}
# Parameter der Verteilung
get_par<- function(vchar,input){
return(
switch(vchar,
"Weibullverteilung"        = { req(input$shape, input$scale); c(input$shape, input$scale) },
"Exponentialverteilung"    = { req(input$rate);               c(input$rate) },
"Gammaverteilung"          = { req(input$shape, input$scale); c(input$shape, input$scale) },
"LogNormalverteilung"      = { req(input$meanlog, input$sdlog); c(input$meanlog, input$sdlog) },
"Weibullverteilung 3 Par." = { req(input$shape, input$scale, input$gamma); c(input$shape, input$scale, input$gamma) }
))
}
# Erneuerungsfunktion
Ht<- function(t, nr_vert, par){
h<- get_moments(nr_vert, par)  #Erwartungswert, Varianz
my<- h[1]
return(t/my+(h[2]-my^2)/(2*my^2))
}
Vt<- function(t, nr_vert, par){
h<- get_moments(nr_vert, par)  #Erwartungswert, Varianz
return(h[2]*t/(h[1]^3))
}
round_up_sig <- function(x) {
magnitude <- 10^floor(log10(x))
ceiling(x / magnitude) * magnitude
}
###################################
dichte <- function(x, nr_vert, par) {
switch(nr_vert,
"1" = dweibull(x, shape = par[1], scale = par[2]),
"2" = dexp(x,    rate  = par[1]),
"3" = dgamma(x,  shape = par[1], scale = par[2]),
"4" = dlnorm(x,  meanlog = par[1], sdlog = par[2]),
"5" = dweibull(x-par[3], shape = par[1], scale = par[2])
)
}
vertfkt <- function(x, nr_vert, par) {
switch(nr_vert,
"1" = pweibull(x, shape = par[1], scale = par[2]),
"2" = pexp(x,    rate  = par[1]),
"3" = pgamma(x,  shape = par[1], scale = par[2]),
"4" = plnorm(x,  meanlog = par[1], sdlog = par[2]),
"5" = pweibull(x-par[3], shape = par[1], scale = par[2])
)
}
# Zuverlaessigkeitsfkt.
zuverlfkt <- function(t, nr_vert, par,tau) {
ii<- ceiling(t/tau)-1     #Index des Intervalls, beginnend mit 0
if (ii==1){
return(1-vertfkt(t, nr_vert, par))
}else{
return(((1-vertfkt(t, nr_vert, par))^ii)*(1-vertfkt(t-ii*tau, nr_vert, par)))
}
}
#Kostenfunktion
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)
}
# tau0_max berechnen
tau0_max<- function(nr_vert, par){
return(fakmaxtau*switch(nr_vert,
"1" = par[2]*gamma(1 + 1/par[1]),
"2" = 1/par[1],
"3" = par[1]*par[2],
"4" = exp(par[1]+par[2]^2/2),
"5" = par[3]+par[2]*gamma(1 + 1/par[1])
))
}
#optimale Präventionszeit berechnen
solve_tau0 <- function(nr_vert, par, interv, c_A, c_P) {
# Linke Seite - Rechte Seite = 0
gleichung <- function(tau0) {
# Integral von (1 - F(x)) dx von 0 bis tau0
integral <- integrate(
function(x){ 1 - vertfkt(x, nr_vert, par)},
lower = 0,
upper = tau0
)$value
h_tau0  <- dichte(tau0, nr_vert, par) /
(1 - vertfkt(tau0, nr_vert, par))  # Hazardrate h(τ₀)
h_tau0 * integral - vertfkt(tau0, nr_vert, par) - c_P / (c_A - c_P)
}
# Nullstelle suchen
uniroot(gleichung, interval = interv,
extendInt = "yes")$root
}
nv  <- get_nr_vert("Weibull")
nv
nv  <- get_nr_vert("Weibullverteilung")
par <- get_par("Weibullverteilung", input) # ggf. manuell: c(2, 100)
par <- c(2, 100)
t_seq <- seq(0, 200, length.out = 500)
S_t <- zuverlfkt(t_seq, nv, par, tau = 0)
for (t in t_seq)zuverlfkt(t, nv, par, tau = 10)
zuverlfkt(t_seq[1], nv, par, tau = 10)
zuverlfkt(t_seq[10], nv, par, tau = 10)
t_seq
t_seq[100]
zuverlfkt(t_seq[100], nv, par, tau = 10)
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')
zuverlfkt <- function(t, nr_vert, par,tau) {
ii<- ceiling(t/tau)-1     #Index des Intervalls, beginnend mit 0
if (ii==1){
return(1-vertfkt(t, nr_vert, par))
}else{
return(((1-vertfkt(tau, nr_vert, par))^ii)*(1-vertfkt(t-ii*tau, nr_vert, par)))
}
}
nr_vert<- 1
par<- c(2,100)
t<- 50
tau<- 30
ii<- ceiling(t/tau)-1
ii
1-vertfkt(t, nr_vert, par)
((1-vertfkt(tau, nr_vert, par))^ii)*(1-vertfkt(t-ii*tau, nr_vert, par))
zuverlfkt(t, nr_vert, par,tau)
runApp('~/Eckhard/stat/R/zuverl/neu')
zuverlfkt <- function(t, nr_vert, par,tau) {
ii<- ceiling(t/tau)-1     #Index des Intervalls, beginnend mit 0
if (ii<=1){
return(1-vertfkt(t, nr_vert, par))
}else{
return(((1-vertfkt(tau, nr_vert, par))^ii)*(1-vertfkt(t-ii*tau, nr_vert, par)))
}
}
runApp('~/Eckhard/stat/R/zuverl/neu')
shiny::runApp("pfad/zur/app")
getwd()
shiny::runApp(getwd())
runApp('~/Eckhard/stat/R/zuverl/neu')
tau <- 50
# Linksseitiger Grenzwert (kurz vor tau)
zuverlfkt(tau - 0.001, nv, par, tau)
# Rechtsseitiger Grenzwert (kurz nach tau)
zuverlfkt(tau + 0.001, nv, par, tau)
# Bei tau selbst
zuverlfkt(tau, nv, par, tau)
1-vertfkt(tau + 0.001, nv, par, tau)
1-vertfkt(tau + 0.001, nv, par)
1-vertfkt(2*tau + 0.001, nv, par)
zuverlfkt(2*tau + 0.001, nv, par, tau)
runApp('~/Eckhard/stat/R/zuverl/neu')
fakmaxprev<- 12 # Faktor fuer praeventive Erneuerungen
runApp('~/Eckhard/stat/R/zuverl/neu')
shiny::runApp(getwd())
shiny::runApp(getwd())
tau <- 50  # dein aktueller Testwert
nv  <- 1
par <- c(2,100)
# Grenzwerte an den Intervallgrenzen
zuverlfkt(0,         nv, par, tau)  # sollte 1 sein
zuverlfkt(tau*0.99, nv, par, tau)  # kurz vor tau
zuverlfkt(tau,       nv, par, tau)  # genau tau
zuverlfkt(tau*1.01,  nv, par, tau)  # kurz nach tau
zuverlfkt(tau*2,     nv, par, tau)  # genau 2*tau
zuverlfkt(tau*2.01,  nv, par, tau)  # kurz nach 2*tau
ceiling(2*tau/tau)-1
ceiling(2.01*tau/tau)-1
zuverlfkt <- function(t, nr_vert, par,tau) {
ii<- ceiling(t/tau)-1     #Index des Intervalls, beginnend mit 0
if (ii<=0){
return(1-vertfkt(t, nr_vert, par))
}else{
return(((1-vertfkt(tau, nr_vert, par))^ii)*(1-vertfkt(t-ii*tau, nr_vert, par)))
}
}
tau <- 50  # dein aktueller Testwert
nv  <- 1
par <- c(2,100)
# Grenzwerte an den Intervallgrenzen
zuverlfkt(0,         nv, par, tau)  # sollte 1 sein
zuverlfkt(tau*0.99, nv, par, tau)  # kurz vor tau
zuverlfkt(tau,       nv, par, tau)  # genau tau
zuverlfkt(tau*1.01,  nv, par, tau)  # kurz nach tau
zuverlfkt(tau*2,     nv, par, tau)  # genau 2*tau
zuverlfkt(tau*2.01,  nv, par, tau)  # kurz nach 2*tau
runApp('~/Eckhard/stat/R/zuverl/neu')
runApp('~/Eckhard/stat/R/zuverl/neu')
runApp('~/Eckhard/stat/R/zuverl/neu')
# Zuverlaessigkeitsfkt.
zuverlfkt <- function(t, nr_vert, par,tau) {
ii<- ceiling(t/tau)-1     #Index des Intervalls, beginnend mit 0
if (ii<=0){
return(1-vertfkt(t, nr_vert, par))
}else{
return(((1-vertfkt(tau, nr_vert, par))^ii)*(1-vertfkt(t-ii*tau, nr_vert, par)))
}
}
tau <- 50  # dein aktueller Testwert
nv  <- 1
par <- c(2,100)
# Grenzwerte an den Intervallgrenzen
zuverlfkt(0,         nv, par, tau)  # sollte 1 sein
zuverlfkt(tau*0.99, nv, par, tau)  # kurz vor tau
zuverlfkt(tau,       nv, par, tau)  # genau tau
zuverlfkt(tau*1.01,  nv, par, tau)  # kurz nach tau
zuverlfkt(tau*2,     nv, par, tau)  # genau 2*tau
zuverlfkt(tau*2.01,  nv, par, tau)  # kurz nach 2*tau
1-vertfkt(tau*2.01,  nv, par, tau)
1-vertfkt(tau*2.01,  nv, par)
runApp('~/Eckhard/stat/R/zuverl/neu')
runApp('~/Eckhard/stat/R/zuverl/neu')
runde_step <- function(x) {
magnitude <- 10^floor(log10(x))
round(x / magnitude) * magnitude
}
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')
max(c(F,F,F,F))
runApp('~/Eckhard/stat/R/zuverl/neu')
fakmaxkosten<- 20
runApp('~/Eckhard/stat/R/zuverl/neu')
nv<- 1
par<- c(2,100)
tau0_max<- tau0_max(nv,par)
hbeg<- step<- tau0max/300
tau0max<- tau0_max(nv,par)
tau0_max<- function(nr_vert, par){
return(fakmaxtau*switch(nr_vert,
"1" = par[2]*gamma(1 + 1/par[1]),
"2" = 1/par[1],
"3" = par[1]*par[2],
"4" = exp(par[1]+par[2]^2/2),
"5" = par[3]+par[2]*gamma(1 + 1/par[1])
))
}
tau0max<- tau0_max(nv,par)
hbeg<- step<- tau0max/300
tau_seq <- seq(hbeg, tau0max, length.out = 300)
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))
j<- 1
kend<- kosten_seq[300]     #Endwert Kosten
pos <- max(which(kosten_seq > max(fakmaxkosten*min(kosten_seq),2*kend))) #Position zum Abschneiden
hbeg<- tau_seq[pos+1]
kosten_seq
max(fakmaxkosten*min(kosten_seq),2*kend)
min(kosten_seq)
runApp('~/Eckhard/stat/R/zuverl/neu')
shiny::runApp(getwd())
runApp('~/Eckhard/stat/R/zuverl/neu')
