######################################
# KB auf R Funktion
# Diese Funtion bestimmt die kritische Konstante q wie in Kapitel 1.3 angegeben.
# Dazu braucht sie die Daten data, um die Anzahl an Beobachtungen nobs zu bestimmen. Man muss ihr ansonsten
# noch das Konfidenznivea alpha und den Grad des Polynoms grad übergeben. Da zur Bestimmung der Werte
# lower und upper inv.X, beta und sigma gebraucht werden, übergebe ich diese auch.
# Rückgabewert sind points,lower,upper,factor. Dabei ist factor der kritische Wert
KB.R <- function(alpha, nobs, grad, inv.X){
# kritischen Wert bestimmen
q <-qf(1-alpha,grad+1, nobs-grad-1)
factor <- sqrt((grad+1)*q)
erg=list(factor)
erg
}
######################################
# KB auf R Funktion pr?f funktion
# Diese Funtion bestimmt die kritische Konstante q, um zu pr?fen ob zwei Modelle gleich sind.
# Dazu brauch sie die Daten data, um die Anzahl an Beobachtungen nobs bestimmen. Man muss ihr ansonsten
# noch das Konfidenznivea alpha und den Grad des Polynoms grad ?bergeben. Da zur Bestimmung der Werte
# lower und upper inv.X, beta und sigma gebraucht werden, ?bergebe ich diese auch.
# R?ckgabewert sind points,lower,upper,factor. Dabei ist factor der kritische Wert
KB.R.pruef <- function(alpha, nobs, grad, k){
# kritischen Wert bestimmen
q <-qf(1-alpha,k,nobs-grad-1)
factor <- sqrt(k*q)
erg=list(factor)
erg
}
##############################
# KB auf Rechteck Funktion
# Beschreibung nach dem Muster:
# Diese Funktion berechnet ein Konfidenzband auf (0,1) und ist dementsprechend nur f?r normierte Daten
# sinnvoll. Eingabewerte sind das Konfidenznivea alpha, die Daten data, die Anzahl an unabh?nigigen
# Parametern grad (Im Hinterkopf ist die Anwendung auf Polynome), die Anzahl an Iterationen niter
# das inverse der Designmatrix inv.X und Sch?tzer f?r beta und sigma
# Die Funktion sollte genau wie in Kapitel 1.4 funktionieren
# Die R?ckgabewerte sind points,lower,upper,factor wie in der letzten Funktion
# Konfidenzband auf (min(x.1), max(x.1)) bestimmen
KB.minmax <- function(alpha, nobs, grad, niter, inv.X, a, b){
# Simulation um den kritischen Wert c zu bestimmen
# niter sollte mindestens 50.000 sein
# feste Werte erzeugen
v <- nobs-grad-1
D <- diag(grad+1)
# Wurzel P aus inv.X=(X'X)^{-1} berechnen
# dazu mache ich eine eigenwertzerlegung von inv.X=ev %*% ew %*% ev.1.
# Dann ist P = ev %*% sqrt(ew) %*% ev.1
eig <- eigen(inv.X)
ew <- diag(length(eig$values))
for(i in 1:length(eig$values)){ew[i,i]=eig$val[i]}
ew.sq <- sqrt(ew)
ev <- eig$vec
ev.1 <- solve(ev)
#ev %*% ew %*% ev.1 # sollte gleich inv.X sein
P <- ev %*% ew.sq %*% ev.1
P.1 <- solve(P)
# Matrix A bestimmen
A.temp=matrix(numeric(),ncol=grad+1)
if(grad>1){
for(i in 1:(grad-1)){A.rows=rbind((D[,i+1]-b*D[,i]) %*% P.1,(a*D[,i]-D[,i+1]) %*% P.1)
A.temp=rbind(A.temp,A.rows)}
A.last <- -D[,1] %*% P.1
A <- rbind(A.temp,A.last)
} else {
A = -D[,1] %*% P.1
}
# Simulation beginnt
temp <- 1
S <- rep(NA, niter)
while(temp<(niter+1))
{
#rmvnorm(n, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)),
# method=c("eigen", "svd", "chol"), pre0.9_9994 = FALSE)
N <- mvtnorm::rmvnorm(1,mean=rep(0,grad+1),D)
sig<-sqrt(rchisq(1,v)/v)
T.vec <- N/sig
#solve.QP(Dmat, dvec, Amat, bvec, meq=0, factorized=FALSE)
sol <- quadprog::solve.QP(Dmat=D,dvec=T.vec,Amat=t(A))
p.T <- norm(sol$solution,type="2")
sol.minus <- quadprog::solve.QP(Dmat=D,dvec=-T.vec,Amat=t(A))
p.T.minus <- norm(sol.minus$solution,type="2")
S[temp]=max(p.T,p.T.minus)
temp <- temp+1
}
# kritische Wert c bestimmen
factor <- quantile(S,1-alpha)
erg=list(factor)
erg
}
##############################
# KB auf Rechteck f?r Polynome im homoskladistischen Fall
#' @useDynLib KBminmaxpoly
#' @importFrom Rcpp sourceCpp
# Berechnet den kritischen Wert des Konfidenzbandes im homoskladistischen Fall. Eingabewerte sind wie oben,
# aber dismal ist niter der grad des Polynoms und nicht nur die Anzahl an unabh?ngigen Parametern.
# Die Funktion funktioniert wie in Kapitel 1.5 nur die Berechnung von den Maxima ist besonders. Daf?r
# sehe die Kommentare im Programmcode.
# R?ckgabewerte wie immer.
KB.poly <- function(alpha, nobs, grad, niter, inv.X, a, b){
# Diese Funktionen braucht man sp?ter in der Berechnung, um das Maximum zu bestimmen
# x.fun berechnet den Wert an x und x.fun.prime an der Ableitung, die f?r Polynome leicht zu
# bestimmen ist
# grad+1, da wir grad des polynoms + konstante brauchen
g.fun.prime <- function (x, T.vec, grad) {(t(x_fun_prime_cpp(x, grad)) %*% t(T.vec)) * (t(x_fun_cpp(x, grad)) %*% inv.X %*% x_fun_cpp(x, grad)) - t(x_fun_cpp(x, grad))%*%t(T.vec) * t(x_fun_prime_cpp(x, grad))%*%inv.X%*%x_fun_cpp(x, grad)}
S.fun = function(v, ngridpoly, xseq, mod, ss){
# Wert simulieren
N <- mvrnormArma(1, rep(0, grad+1), inv.X)
sigma.hat <- sqrt(rchisq(n=1,df=v)/v)
# Berechnug der Maxima
# orientiert sich an der funktion uniroot.all aus dem rootSolve Package
# uniroot.all does that by first subdividing the interval into small sections and, for all sections
# where the function value changes sign, invoking uniroot to locate the root.
#S_fun_cpp(x, grad, sigma.hat, N, inv.X)
g.lapply=function(x){g.fun.prime(x, T.vec, grad)}
mod=sapply(xseq, g.lapply)
# ist mod==0 liegt schonmal eine Wurzel vor
roots.1 <- xseq[which(mod == 0)]
# multipliziert aufeinanderfolgende Funktionswerte
ss <- mod[1:ngridpoly] * mod[2:(ngridpoly + 1)]
# ist ss[i]<0, so liegt ein Vorzeichenwechsel vor und zwischen mod[i] und mod[i+1]
# liegt ein Nullstelle -> ii speichert die Inidizes der Werte von ss, die kleiner Null sind
ii <- which(ss < 0)
# loop ?ber alle indizes in ii. Bei jedem Durchlauf Aufruf der uniroot-Funktion, die Wurzeln im
# Eindimensionalen findet. Alle diese Wurzeln werden in all.roots gespeichert.
# wrap-around function, da uniroot nur eindimensional funktionen vertr?gt
roots.2=rep(NA, length(ii))
temp.1=1
for (i in ii){roots.2[temp.1]= (uniroot(function(x)(g.fun.prime(x,T.vec, grad)), lower = xseq[i], upper = xseq[i + 1])$root)
temp.1=temp.1+1}
all.roots=c(roots.1,roots.2)
# als n?chstes soll g.fun(i) f?r i in all.roots gefunden werden
# Fallunterscheidung, ob Wurzeln gefunden wurden
if(is.na(all.roots[1])==T){erg=max(abs(g.fun(a, T.vec, grad)),abs(g.fun(b, T.vec, grad)))}
else {
# all.g speichert die Funktionswerte von g.fun an den Nullstellen von g.fun.prime
all.g <- rep(NA, length(all.roots))
temp=c(1:length(all.roots))
g.fun.apply=function(x){g_fun_cpp(x, grad, T.vec, inv.X)}
all.g=sapply(temp, g.fun.apply)
erg=max(abs(g_fun_cpp(a, grad, T.vec, inv.X)),abs(g_fun_cpp(b, grad, T.vec, inv.X)),abs(all.g))
}
erg
}
# feste Werte erzeugen
v=nobs-grad-1
ngridpoly = 100 # Feinheit des Grids festlegen, auf dem wir nachher die wurzeln von g.fun.prime suchen
xseq <- seq(a, b, len = ngridpoly + 1) # subdividing the interval [a,b] in ngrid+1 parts
mod = rep(NA,(ngridpoly+1)) # mod speichert die Funktionswerte auf [a,b] mit Abstand ngrid
ss=rep(NA,(ngridpoly+1)) # multipliziert aufeinander folgende Funktionswerte
S <- rep(NA, niter) # lehrer Vektor f?r die Simulationsergebnisse
# Simulation durchf?heren
for(i in 1:niter)
{
S[i]=S.fun(v, ngridpoly, xseq, mod, ss)
}
# kritischen Wert bestimmen
factor <- quantile(S,1-alpha)
erg=list(factor)
erg
}
#################################################
# Diese Funktion bestimmt das Maximum der Funktion g auf einem Grid mit Feinheit ngridpoly auf [0,1]
# indem die Funktionswerte von jedem Wert bestimmt werden.
# Dies müsste schneller, als der letzte Ansatz funktionieren
KB.poly.fast <- function(alpha, nobs, grad, niter, inv.X, a, b, ngridpoly){
# von dieser Funktion m?chte man das Maximum bestimmen
S.fun = function(v, values, xseq, i){
# Wert simulieren
sigma.hat = sqrt(rchisq(n=1,df=v)/v)
N = matrix(mvtnorm::rmvnorm(1,mean=rep(0,grad+1),inv.X), ncol=grad+1)
# Berechnug der Maxima
# Bestimmt einfach den Wert von g.fun auf einem Grid auf [0,1] mit Feinheit ngridpoly
g.lapply=function(x){S_fun_cpp(x, grad, sigma.hat, N, inv.X)}
values=rep(NA, length(xseq))
values=sapply(xseq, g.lapply)
# Erzeugt plots der S_fun_cpp Funktion. In der Ausarbeitung wird diese Funktion mit h bezeichnet
#pdf(file = paste("tests/Funktion_h/test_",i,".pdf", sep=""))
#plot(xseq,values)
#dev.off()
# Maximum der Funktion g in dieser Simulation
erg=max(abs(values))
erg
}
# feste Werte erzeugen
v=nobs-grad-1
xseq <- seq(a, b, len = ngridpoly + 1) # subdividing the interval [a,b] in ngrid+1 parts
values = rep(NA,(ngridpoly+1)) # mod speichert die Funktionswerte auf [a,b] mit Abstand ngrid
S <- rep(NA, niter) # lehrer Vektor f?r die Simulationsergebnisse
# Simulation durchf?heren
for(i in 1:niter){
S[i]=S.fun(v,values, xseq, i)
}
# kritischen Wert bestimmen
factor <- quantile(S,1-alpha)
erg=list(factor)
erg
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.