all.reg <- function (X,y) 
{## X on regressioanalyysin X-matriisi ilman ykkössaraketta
 ## y on vaste

 n <- nrow(X)
 bn <- log(n)
 p <- ncol(X)
 k <- 2^p-1
 info <- matrix(0,k+1,5)
 colnames(info) <- c("AIC","BIC","Cp","PRESS","No of var's")
 g <- lm(y~1)
 gf <- lm(y~X)
 RSS <- sum(g$res^2)
 sigma2 <- sum(gf$res^2)/(n-p-1) 
 info[1,1:2] <- n*log(RSS/n)
 info[1,3] <- RSS/sigma2 - n
 info[1,4] <- sum((g$res/(1-influence(g)$hat))^2)
 info[1,5] <- 0
 for (i in 2:(k+1)) { b <- binary(i-1)
     g <- lm(y~X[,c(b$dico,rep(F,p))[1:p]])
     RSS <- sum(g$res^2)
     info[i,1:2] <- n*log(RSS/n) + sum(b$binary)*c(2,bn)
     info[i,3] <- RSS/sigma2 + 2*sum(b$binary) - n
     info[i,4] <- sum((g$res/(1-influence(g)$hat))^2)
     info[i,5] <- sum(b$binary)
                 }
  out <- list(info=info, X=X, y=y)
  out
}

##################################################################

best.reg <- function (info) 
{## info on output funktiosta all.reg

 p <- ncol(info$X)
 best <- array(" ",c(4,p))
 colnames(best) <- colnames(info$X)
 rownames(best) <- colnames(info$info)[1:4]
 n <- rep(0,4)
n[1]<-which.min(info$info[,"AIC"])
n[2] <-which.min(info$info[,"BIC"])
n[3]<-which.min(info$info[,"Cp"])
n[4]<-which.min(info$info[,"PRESS"])
 best[1,c(binary(-1+n[1])$dico,rep(F,p))[1:p]] <- "*"
 best[2,c(binary(-1+n[2])$dico,rep(F,p))[1:p]] <- "*"
 best[3,c(binary(-1+n[3])$dico,rep(F,p))[1:p]] <- "*"
 best[4,c(binary(-1+n[4])$dico,rep(F,p))[1:p]] <- "*"

 print(best, quote=F)
 print(info$info[n,])
 invisible(best)
}

#####################################################################

outlier.test <- function (lm.obj, level=0.05) 
{## Finds outliers using Bonferroni's and Hochberg's method, 
 ## Biometrika 1988, 800-2
 ## lm.obj on output funktiosta lm ja level on merkitsevyystaso

 res <- lm.obj$resid
 n <- length(res)
 df <- lm.obj$df
 sigma <- summary(lm.obj)$sigma
 leverage <- 1 - influence(lm.obj)$hat
 r <- res/(sigma*sqrt(leverage))
 tres <- r*sqrt((df-1)/(df-r^2))
 pval <- 2*(1-pt(abs(tres),df=df))
#print(round(pval,4))
 print("Outliers found:")
 print("Bonferroni")
 b <- (1:n)[pval <= level/n]
 if (is.null(b)) print("None") else 
       {outlier <- cbind(tres[b],pval[b])
        colnames(outlier) <- c("t-arvo","Bonf. p-arvo")
        print(outlier)
       }

 print("Hochberg")
 h <- (1:n)[pval <= level/(n + 1 - rank(pval))]
 if (is.null(h)) print("None")
    else {m <- max(pval[h])
          h <- (1:n)[pval <= m]
          print(h)
         }
# print(round(cbind(pval,level/(n + 1 - rank(pval))),4))

}

######################################################################
diagtest <- function(lm.out,z=NULL,test="Linearity",repl=10){
 ## Possible tests:
 ## "Linearity", "HetSced", "Autocorrelation", "Normality"
 ##  
 ## lm.out on lm- objekti, so. lm-funktion tulos
 ## z = alternative variables for "Linearity" and "HetSced"

Imhof <- function (x,ev,df=1,ncp=NULL) 
{## Computes probabilities P(sum(ev[j]*Chi_j^2) > x) using Imhof's   technique 
 ## where Chi_j^2 's are independent Chi-squared variables with df[j] 
 ## degrees of freedom and with noncentrality parameters np[j].
 ## Default df's are = 1, and noncentrality parameters = 0

 s<-mean(abs(ev))
 x<-x/s
 ev<-ev/s
 
 CF <- function(u,x=x,ev=ev,df=df,ncp=ncp) ##
  {uev<-outer(u,ev,"*")
   uev2<-uev*uev
   if(length(df)==1) {theta <-apply(atan(uev),1,"sum")*df
                      rho<-(apply(1+uev2,1,"prod"))^(df*.25)}
   
   else {theta <- atan(uev)%*%df
         power <- t(array(0.25*df,c(length(df),length(u))))
         rho  <- apply((1+uev2)^power,1,"prod")}

 if(!is.null(ncp)) {theta <-theta+(uev/(1+uev2))%*%ncp
                   rho<-rho*exp(.5*(uev2/(1+uev2))%*%ncp)}
 
 theta <-sin(.5*(theta-x*u))
 ratio <-theta/rho
 ratio <-(u!=0)*ratio+(u==0)*0.5*(sum(ev*(df+ncp))-x)
 ratio <-c(ratio/(u+(u==0)))
 return(ratio)
 }
 integ <- integrate(CF,0,Inf,x=x,ev=ev,df=df,ncp=ncp)
 prob <- 0.5 + integ$value/pi
 out <- c(list(prob=prob),integ)
 return(out)
 }

 X <- model.matrix(lm.out)
 p <- ncol(X)
 n <- nrow(X)
 res <- lm.out$res
 Warn <- " "
 names(Warn) <- " "
 Text2 <- Warn


 if (test=="Linearity") {Text <- "Test for linearity:"
 W <- outer(1:n,1:n,"pmin")
 
 L <- pval <- ave <- evar <- gappr <- NULL

 if (is.null(z)) z <- as.matrix(X[,-1])
 q <- ncol(as.matrix(z))
 for (j in 1:q) {r <- order(as.vector(z[,j]))
 zr <- as.vector(z[r,j])
 Xr <- X[r,]; resr <- res[r]
 W <- outer(zr,zr,"pmin")
 Q <- diag(n)-Xr%*%solve(t(Xr)%*%Xr)%*%t(Xr)
 L <- c(L,c(resr%*%W%*%resr)/((zr[n]-zr[1])*sum(resr*resr)))
 eig <- eigen(Q%*%W%*%Q,symmetric=T,only.values=T)
 pval <- c(pval,Imhof(x=0,ev=eig$val[1:(n-p)]/(zr[n]-zr[1]) - L[j])$prob)

 ave <- c(ave,mean(eig$val[1:(n-p)]/n))
 evar <- c(evar,2*mean((eig$val[1:(n-p)]/n - ave[j])^2/(n-p+2)))
 gappr <- c(gappr,1-pgamma(L[j],shape=ave[j]^2/evar[j],scale=evar[j]/ave[j]))
 }

 out <- cbind(L,pval)
 colnames(out) <- c("Test value","Pr(>Test value)")
 rownames(out) <- colnames(X)[2:p]

 if (q > 1) {r <- order(lm.out$fitted)
 Xr <- X[r,]; resr <- res[r]
 zr <- sort(lm.out$fitted)
 W <- outer(zr,zr,"pmin")
 Q <- diag(n)-Xr%*%solve(t(Xr)%*%Xr)%*%t(Xr)
 L <- c(L,c(resr%*%W%*%resr)/((zr[n]-zr[1])*sum(resr*resr)))
 eig <- eigen(Q%*%W%*%Q,symmetric=T,only.values=T)
 pval <- c(pval,Imhof(x=0,ev=eig$val[1:(n-p)]/(zr[n]-zr[1]) - L[q+1])$prob)

 out <- cbind(L,pval)
 colnames(out) <- c("Test value","Pr(>Test value)")
 rownames(out) <- c(colnames(X)[2:p],"Fitted values")
            }
  }
 
 if (test=="HetSced") {Text <- "Test for constant variance:"
 
 L <- pval <- ave <- evar <- gappr <- NULL

 if (is.null(z)) z <- as.matrix(X[,-1])
 q <- ncol(as.matrix(z))
 for (j in 1:q) {W <- as.vector(z[,j]); W <- W - min(W) +1 
 Q <- diag(n)-X%*%solve(t(X)%*%X)%*%t(X)
 L <- c(L,sum(W*res^2)/(n*sum(res*res))); k <- length(L)
 eig <- eigen(Q%*%diag(W)%*%Q,symmetric=T,only.values=T)
 pval <- c(pval,Imhof(x=0,ev=eig$val[1:(n-p)]/n - L[k])$prob)
 ave <- c(ave,mean(eig$val[1:(n-p)]/n))
 evar <- c(evar,2*mean((eig$val[1:(n-p)]/n - ave[k])^2/(n-p+2)))
 gappr <- c(gappr,1-pgamma(L[k],shape=ave[k]^2/evar[k],scale=evar[k]/ave[k]))
 }

 out <- cbind(L,pval,1-pval,2*pmin(pval,1-pval))
 colnames(out) <- c("Test value","Pr(>Test value)","Pr(<Test value)",
                    "Two-sided")
 rownames(out) <- colnames(X)[2:p]
 
 if (q > 1){
 Q <- diag(n)-X%*%solve(t(X)%*%X)%*%t(X)
 W <- lm.out$fitted
 W <- W - min(W) +1
 L <- c(L,sum(W*res^2)/(n*sum(res*res)))
 eig <- eigen(Q%*%diag(W)%*%Q,symmetric=T,only.values=T)
 pval <- c(pval,Imhof(x=0,ev=eig$val[1:(n-p)]/n - L[length(L)])$prob)
 out <- cbind(L,pval,1-pval,2*pmin(pval,1-pval))
 colnames(out) <- c("Test value","Pr(>Test value)","Pr(<Test value)",
                    "Two-sided")
 rownames(out) <- c(colnames(X)[2:p],"Fitted values")
             }

 Text2 <- rbind(c("2nd col. = p value against increasing variance"), 
                c("3rd col. = p value against decreasing variance"),
                c("4th col. = p value against increasing or decreasing variance"))
 colnames(Text2) <- " "
 rownames(Text2) <- c(" "," "," ")
 Warn <- "CAUTION: Normality assumption crucial!!!"
 names(Warn) <- " "
 }

 if (test=="Autocorrelation") {Text <- "Durbin-Watson test for no autocorrelation:"
 W <- 2*diag(n) - (abs(outer(1:n,1:n,"-"))==1)
 W[1,1] <- W[n,n] <- 1
 L <- sum(diff(res)^2)/sum(res^2)
 Q <- diag(n)-X%*%solve(t(X)%*%X)%*%t(X)
 eig <- eigen(Q%*%W%*%Q,symmetric=T,only.values=T)
 pval <- 1-Imhof(x=0,ev=eig$val[1:(n-p)] - L)$prob
 ave <- mean(eig$val[1:(n-p)])
 evar <- 2*mean((eig$val[1:(n-p)] - ave)^2/(n-p+2))
 gappr <- pgamma(L,shape=ave^2/evar,scale=evar/ave)

 out <- cbind(L,pval,1-pval,2*pmin(pval,1-pval))
 colnames(out) <- c("Test value","Pr(<Test value)","Pr(>Test value)",
                    "Two-sided")
 rownames(out) <- " "
 Text2 <- rbind(c("2nd col. = p value against positive autocorrelation"), 
                c("3rd col. = p value against negative autocorrelation"),
                c("4th col. = p value against both positive and negative autocorrelation"))
 colnames(Text2) <- " "
 rownames(Text2) <- c(" "," "," ")
 }

 if (test=="Normality") {Text <- "Anderson-Darling test for normality via simulation:"
 res <- sort(res); res <- res/sqrt(sum(res*res)/(n-p))
 w <- 2*(1:n)-1
 pr <- pnorm(res)
 AD <- -mean(w*log(pr*(1-pr[n:1]))) - n
 ADs <- NULL
 for (i in 1:repl) {
 y <- rnorm(n)
 ou <- lm(y~X-1)
 ress <- sort(ou$res)
 ress <- ress/sqrt(sum(ress*ress)/(n-p))
 prs <- pnorm(ress)
 ADs <- c(ADs,-mean(w*log(prs*(1-prs[n:1]))) - n)
     }
 pval <- mean(ADs>AD)
 bound <- qnorm(.975)*sqrt(pval*(1-pval)/repl)
 out <- cbind(AD,pval,pval-bound,pval+bound,repl)
 rownames(out) <- " "
 colnames(out) <- c("Test value","Pr(>Test value)","2.5%","97.5%",
                    "No of repl's")
 }

 Call <- "Call:"
 names(Text) <- " "
 names(Call) <- " "
 
 out <- list(Call=Call,call=lm.out$call,Text=Text,Table=out,Text2=Text2,
             Warn=Warn)

 invisible(out)
 
 print(out$Call,quote=F)
 print(out$call)
 print(out$Text,quote=F)
 print(out$Table,print.gap=3)
 print(out$Text2,quote=F)
 print(out$Warn,quote=F)
}

#########################################################################
exp.reg<-function (y,X,par) 
{# Estimoi ILS-menetelmällä eksponenttisen regression parametrit

g <- function(X,par){
   g0 <- par[1]+par[2]*exp(par[3]*X)
   g0
                     }
 g.der <- function(X,par){
    Z0 <- cbind(1,exp(par[3]*X),par[2]*X*exp(par[3]*X))
    Z0
                          }
 err <- 1
 iter <- c(par,sum((y-g(X,par))^2))
 names(iter) <- c("beta","gamma","delta","RSS")
 while(err > 1e-6) {
  g0 <- g(X,par)
  Z0 <- g.der(X,par)
  delta <- qr.solve(Z0,y-g0)
  par <- par+delta
  err <- max(abs(delta))
  iter <- rbind(iter,c(par,sum((y-g(X,par))^2)))
                   }
 print(iter)
 out <- par
}