library(MASS)

coord=function(xlim=c(0,5),ylim=c(0,5)){
   x=0
   y=0
   plot(x,y,xlim=xlim,ylim=ylim,type="n")
   }
   
shift=function(a,from=c(0,0),text=deparse(substitute(a))){
   arrows(from[1],from[2],from[1]+a[1],from[2]+a[2])
   text(locator(1),text)
   }
   
point=function(a,text=deparse(substitute(a))){
   points(a[1],a[2],cex=2)
   text(locator(1),text)
   }

line=function(from=c(0,0),to=c(0,1),text=NULL,...){
   line(from[1],from[2],to[1],to[2],...)
   text(locator(1),text)
   }

rand.matrix=function(m,n,rank=NULL){
   if (is.null(rank)) rank=min(m,n)
   a=matrix(rnorm(m*rank),m,rank)
   b=matrix(rnorm(rank*n),rank,n)
   return(a%*%b)
   }

unim.matrix=function(m){
   a=matrix(round(3*runif(m*m))-2,m,m)
   diag(a)=rep(1,m)
   b=a
   a[upper.tri(a)]=0
   b[lower.tri(b)]=0
   return(a%*%b)
   }
   
rand.orth=function(d){
   x=randn(d,d)
   return(qr.Q(qr(x)))
   }
   
rand.cov=function(d){
   x=rand.orth(length(d))
   return(x%*%diag(d)%*%t(x))
   }      

rand.stocks=function(n,d,mu=rep(0,d),cov=diag(rep(1,d))){
   mu=mu/360
   cov=cov/360
   x=mvrnorm(n,mu-diag(cov)/2,cov)
   y=apply(x,2,cumsum)
   z=rbind(rep(1,d),exp(y))
   return(ts(z,deltat=1/360,start=0))
   }

im=function(a){
   tmp=qr(a)
   return(qr.Q(tmp)[,1:tmp$rank])
   }
   
ker=function(a){Null(t(a))}   

surface=function(from=0,to=5,frmla,...){
   x=seq(from,to,len=30)
   y=x
   z=outer(x,y,function(x,y) eval(parse(text=frmla)))
   persp(x,y,z,ticktype="detailed",...)
   }

grad=function(x,y,frmla,shr=0.01){ 
   vx=as.vector(matrix(x,length(x),length(y)))
   vy=as.vector(matrix(y,length(x),length(y),byrow=TRUE))
   x=vx
   y=vy
   df=deriv(parse(text=frmla),c("x","y"))
   grad=attr(eval(df),"gradient")
   return(list(x0=vx,y0=vy,x1=vx+shr*grad[,1],y1=vy+shr*grad[,2]))
   }


gfield=function(from=-3,to=3,frmla="x^2+y^2"){
   x=seq(from,to,len=20)
   y=x   
   coord(xlim=c(from,to),ylim=c(from,to))   
   g=grad(x,y,frmla,shr=0.2)
   arrows(g$x0,g$y0,g$x1,g$y1,length=0.05)
   z=outer(x,y,function(x,y) eval(parse(text=frmla)))
   contour(x,y,z,add=TRUE)
   }



norm=function(x){
   return(sqrt(sum(x^2)))
   }

rref=function(x,integer=FALSE){
   m=dim(x)[1]
   n=dim(x)[2]
   y=diag(rep(1,m))
   i=1
   j=1
   repeat {
      while (abs(x[i,j])<1.0e-10) {
          temp=x[,j]
          temp[1:i]=0
          k=order(abs(temp))[m]
          if (abs(temp[k])<1.0e-10) {
             j=j+1
             if (j>n) break
             }
          else {
             x[i,]=x[i,]+x[k,]
             y[i,]=y[i,]+y[k,]
             }
          }   
      if (j>n) break   
      if (integer==FALSE) {
              y[i,]=y[i,]/x[i,j] 
              x[i,]=x[i,]/x[i,j] 
      }
      y[-i,]=x[i,j]*y[-i,,drop=FALSE]-x[-i,j,drop=FALSE]%*%y[i,,drop=FALSE]      
      x[-i,]=x[i,j]*x[-i,,drop=FALSE]-x[-i,j,drop=FALSE]%*%x[i,,drop=FALSE]  
      i=i+1
      j=j+1
      if ((i>m)|(j>n)) break
      }
   if (integer==FALSE){   
   for (i in (1:m)){
      j=which(!(abs(x[i,])<1.0e-10))   
      if (length(j)>0) {
         k=min(j)
         y[i,]=y[i,]/x[i,k]
         x[i,]=x[i,]/x[i,k]
         }
      else {
         j=which(!(abs(y[i,])<1.0e-10))
         if (length(j)>0) {
            k=min(j)
            y[i,]=y[i,]/y[i,k]
            }
         }   
      }}   
   return(list(echelon.form=round(x,8),row.operations=round(y,8)))
   }   

info2=function(y,x1,x2){
   r.yx1=cor(y,x1)
   r.yx2=cor(y,x2)
   r.yx1.x2=cor(lm(y~x2)$residuals,lm(x1~x2)$residuals)
   r.yx2.x1=cor(lm(y~x1)$residuals,lm(x2~x1)$residuals)
   return(list(cor=c(r.yx1,r.yx2),pcor=c(r.yx1.x2,r.yx2.x1),
          rsquare=cbind(c(r.yx1^2,r.yx2^2),
                        c(r.yx2.x1^2*(1-r.yx1^2),r.yx1.x2^2*(1-r.yx2^2)))))
   }                     
   
   
   
# Eckpunkte des d-dimensionalen Einheitswürfels

cube <- function(d){
   x <- 0:(2^d-1)
   cube <- numeric()
   while(any(x>0)){
      temp <- floor(x/2)
      cube <- rbind(cube,x-temp*2)
      x <- temp
      }
   return(t(cube))
   }

# Modellgraph bei multipler Regression

modelselection <-  function(xx){
   n.data <- dim(xx)[1]
   n.pred <- dim(xx)[2]-1
   response <- xx[,1]
   predictors <- xx[,-1]
   models <- cube(n.pred)
   n.models <- dim(models)[1]
   models <- models[n.models:1,]
   out <- numeric()
   for (i in 1:n.models){
      design <- cbind(rep(1,n.data),predictors[,models[i,]==1])
      res <- response-design%*%(ginv(design)%*%response)
      ssr <- sum(res*res)
      if (i==1) {ssr0 <- ssr}
      if (i<n.models) 
         rsquare <- round((1-ssr/(var(response)*(n.data-1))),5)
      else
         rsquare=""    
      n.vars <- sum(models[i,])
      if (i>1)
         {fv <- (ssr-ssr0)/ssr0*(n.data-n.pred-1)/(n.pred-n.vars)
          pv <- round((1-pf(fv,n.pred-n.vars,n.data-n.pred-1)),5)}
         else
         {pv <- ""}   
      out <- rbind(out,c(n.vars,rsquare,pv,models[i,]))
      }
    i <- order(out[,1])
    out <- out[i,]
    out <- out[n.models:1,]  
    out <- data.frame(out)
    names <- attr(xx,"dimnames")[[2]][2:(n.pred+1)]
    attr(out,"names") <- c("ExplVar","RSquare","PValue",names)
    return(out)
    }




modelselection.1 <-  function(xx){
   n.data <- dim(xx)[1]
   n.pred <- dim(xx)[2]-1
   response <- xx[,1]
   predictors <- xx[,-1]
   models <- cube(n.pred)
   n.models <- dim(models)[1]
   models <- models[n.models:1,]
   out <- numeric()
   pts=numeric(0)
   lbs=character(0)
   for (i in 1:n.models){
      design <- cbind(rep(1,n.data),predictors[,models[i,]==1])
      res <- response-design%*%(ginv(design)%*%response)
      ssr <- sum(res*res)
      if (i==1) {ssr0 <- ssr}
      if (i<n.models) 
         rsquare <- round((1-ssr/(var(response)*(n.data-1))),5)
      else
         rsquare=""    
      n.vars <- sum(models[i,])
      if (i>1)
         {fv <- (ssr-ssr0)/ssr0*(n.data-n.pred-1)/(n.pred-n.vars)
          pv <- round((1-pf(fv,n.pred-n.vars,n.data-n.pred-1)),5)}
         else
         {pv <- ""}   
      if ((pv>0.01)|(pv=="")) {   
      if (i<n.models) {
         pts=rbind(pts,c(n.vars,rsquare)) 
         lbs=c(lbs,paste(dimnames(predictors)[[2]][which(models[i,]==1)],collapse=","))
         }
      out <- rbind(out,c(n.vars,rsquare,pv,models[i,]))}
      }
    i <- order(out[,2])
    out <- out[i,]
    out <- out[dim(out)[1]:1,]  
    out <- data.frame(out)
    names <- attr(xx,"dimnames")[[2]][2:n.pred]
    attr(out,"names") <- c("ExplVar","RSquare","PValue",names)
    plot(pts[,1],pts[,2],xlim=c(0,n.pred+1),type="p",pch=19,col="red")
    grid(nx=NA,ny=NULL,col="black")
    identify(pts[,1],pts[,2],lbs)
    return(out)
    }

scatter2=function(x,y,...){
   points(x,y,...)
   abline(coef(lm(y~x)))
   }

rqf=function(type="posdef"){
   repeat {
      x=floor(10*runif(3))-4
      a=x[1]
      d=x[1]*x[3]-x[2]*x[2]
      if ((type=="posdef")&(a>0)&(d>0)) break
      else if ((type=="negdef")&(a<0)&(d>0)) break 
      else if ((type=="indef")&(d<0)) break     
      else if ((type=="semidef")&(d==0)) break   
      }
   return(matrix(c(x[1],x[2],x[2],x[3]),2,2))
   }      

permute <- function(x, byrow=FALSE){
   if (byrow==TRUE) x  <-  t(x)
   m  <-  dim(x)[1]
   n  <-  dim(x)[2]
   temp  <-  as.vector(matrix(1:n,m,n,byrow=TRUE))
   perm  <-  order(runif(m*n)+temp)
   x  <-  matrix(as.vector(x[perm]),m,n)
   if (byrow==TRUE) x  <-  t(x)
   return(x)
   }

rportfolio=function(n,mean=m,cov=c){
   k=length(mean)
   w=runif(n)
   dim(w)=c(n,1)
   for (i in 2:k) w=cbind(w,(1-apply(w,1,sum))*runif(n))
   w=permute(w,byrow=TRUE)
   w=sweep(w,1,apply(w,1,sum),"/")
   mm=w%*%mean
   ss=apply(w*(w%*%cov),1,sum)
   return(list(s=ss,m=mm))
   }