myfunctions.r

## logit/expit
logit = function(x,b=1) log(x/(b-x))
expit = function(x,b=1) b/(1+exp(-x))

## matlab equivalent of repmat in R
repmat = function(X,m,n){
  ## tiles matrix X, m by n times
  mx = dim(X)[1]
  nx = dim(X)[2]
  matrix(t(matrix(X,mx,nx*n)),mx*m,nx*n,byrow=T)
}

acfi = function(x, lag.max = NULL,
         type = c("correlation", "covariance", "partial"),
         plot = TRUE,  demean = TRUE,...)
{
acf(x, lag.max = lag.max,
         type = type,
         plot = plot, na.action = na.pass, demean = TRUE, ...)
}

"%&%" <- function(a, b) paste(a, b, sep="")


rangi = function(x){
 range(x, na.rm = T)
}

cori = function(x,y=NULL){
cor(x, y = y, use = "complete.obs")
}


## myoptim
## minimizes fun using optim with nelder mead
## algorithm. 

myoptim = function(param,fun,epsi=1e-5,hess=F){
  dife = 10
  eini = fun(param)
  while(dife>epsi){
    res=optim(param,fun,hessian=hess)
    dife = abs(eini-res$value)
    param=res$par
    eini=res$value
  }
  res
}

##======================
## mean squared error sum((x-y)^2)/n
mse <- function(x,y) {
  n <- length(x)
  sum((x-y)^2)/n}
                                                                                
rmse <- function(x,y) sqrt(mse(x,y))
##======================
regress <- function(form,plot=F){
  fit=lm(form)
  fitsum=summary(fit);
  print(fitsum$coeff);
  print(c("R squared","residual se"))
  print(round(c(fitsum$r.squared,rmse(fit$res,0)),3))
  }

## Great circle distance from library fields, default is km.
gcdist <- function (loc1, loc2, miles = FALSE, R = NULL) 
{
    if (is.null(R)) {
        if (miles) 
            R <- 3963.34
        else R <- 6378.388
    }
    if (missing(loc2)) 
        loc2 <- loc1
    coslat1 <- cos((loc1[, 2] * pi)/180)
    sinlat1 <- sin((loc1[, 2] * pi)/180)
    coslon1 <- cos((loc1[, 1] * pi)/180)
    sinlon1 <- sin((loc1[, 1] * pi)/180)
    coslat2 <- cos((loc2[, 2] * pi)/180)
    sinlat2 <- sin((loc2[, 2] * pi)/180)
    coslon2 <- cos((loc2[, 1] * pi)/180)
    sinlon2 <- sin((loc2[, 1] * pi)/180)
    pp <- cbind(coslat1 * coslon1, coslat1 * sinlon1, sinlat1) %*% 
        t(cbind(coslat2 * coslon2, coslat2 * sinlon2, sinlat2))
    R * acos(ifelse(pp > 1, 1, pp))
}
## cord distance, default is km.
## loc1 = [lon, lat ]
cordist <- function (loc1, loc2, miles = FALSE, R = NULL) 
{
    if (is.null(R)) {
        if (miles) 
            R <- 3963.34
        else R <- 6378.388
    }
    if (missing(loc2)) 
        loc2 <- loc1
    coslat1 <- cos((loc1[, 2] * pi)/180)
    sinlat1 <- sin((loc1[, 2] * pi)/180)
    coslon1 <- cos((loc1[, 1] * pi)/180)
    sinlon1 <- sin((loc1[, 1] * pi)/180)
    coslat2 <- cos((loc2[, 2] * pi)/180)
    sinlat2 <- sin((loc2[, 2] * pi)/180)
    coslon2 <- cos((loc2[, 1] * pi)/180)
    sinlon2 <- sin((loc2[, 1] * pi)/180)
    pp <- cbind(coslat1 * coslon1, coslat1 * sinlon1, sinlat1) %*% 
        t(cbind(coslat2 * coslon2, coslat2 * sinlon2, sinlat2))
    R * sqrt( 2* (1-ifelse(pp > 1, 1, pp) ) )
}



## matern function with Handcock's parameterization
## hmatern
hmatern = function(u,range,smoo){
  if (is.vector(u)){
    names(u) <- NULL
    y = u
    y[]=1
  }
  if (is.matrix(u)) {
    dimnames(u) <- list(NULL, NULL)
    y = u
    y[,]=1
  }
  unonzero= u>0
  if(smoo<100) {
    urange = u*2*sqrt(smoo)/range
    y[unonzero] = (urange[unonzero]^smoo) * besselK(urange[unonzero],smoo)/
      (2^(smoo - 1))/ gamma(smoo)}
  else
    {y[unonzero] = exp(-u[unonzero]^2)}
  y
  }
##===========================================================
## simple useful functions
##===========================================================
suma <- function(x) {
  indnna = !is.na(x)
  if(sum(indnna)>0) res = sum(x[indnna]) else res=NA
  res
}
maxi <- function(x){
  indnna = !is.na(x)
  if(sum(indnna)>0) res=max(x[indnna]) else res=NA
  res}

mini <- function(x){
  indnna = !is.na(x)
  if(sum(indnna)>0) res=min(x[indnna]) else res=NA
  res}

meana <- function(x){
  indnna = !is.na(x)
  if(sum(indnna)>0) res=mean(x[indnna]) else res=NA
  res}

vari <- function(x){
  indnna = !is.na(x)
  if(sum(indnna)>0) res=var(x[indnna]) else res=NA
  res
}

longi <- function(data){length(data[,1])}

plotfit <- function(fit)
         {
           par( mfrow=c(2,2) )
           plot(fit)
           par( mfrow=c(1,1) )
         }

algunos <- function(x,cant=10) x[1:cant,]

## returns the row and column of an index matrix == 1
indmat <- function(x)
{
  n <- length(x[,1])
  tm1 <- seq(1,n)
  tm2 <- matrix(rep(1,n*n),n)
  rowmat <- diag(tm1)%*%tm2
  colmat <- tm2%*%diag(tm1)
  ro <- rowmat[x==1]
  co <- colmat[x==1]
  cbind(ro,co)
}

reducing margins in R plots

reduce margins
par(mar=c(2,2,1,1)) # default in R is c(5, 4, 4, 2) + 0.1
c(bottom, left, top, right)

confidence intervals ACF

under independence assumption: clim0 = qnorm((1+.95)/2)/sqrt(N)
ARIMA: at k-th lag: clim0 * sqrt( 1+2 \sum_1^k( acf(i)^2 ) ) 

rules of thumbs for ARIMA model selection

 Prof. R Nau's lecture notes

theoretical ACF of ARMA

ARMAacf(ar = numeric(0), ma = numeric(0), lag.max = r, pacf = FALSE)

plot in log scale with axis

plot(1:100, log = "y", yaxt = "n") # do not show y axis
axis(2, c(1,10,100)) # draw y axis with required labels 

Other useful R functions blogs and pages

Paul Johnson's R Page

Li Chen's R page

James Forester R Functions' Blog 

R blog

calculating angle from xy coordinates

## from north (-pi,pi) positive toward east
angle.north = function(xx,yy){sign(xx)*(yy<0)*pi + atan(xx/yy)}

## from east (0,2.pi) positive toward north
angle = function(xx,yy){atan(yy/xx) - sign(xx)*sign(yy)*pi*(xx<0)}

empirical distributions CDF

"empirical CDF's are not good at showing differences in the tails of a distribution"
Owen, A. B., Empirical Likelihood, 2001, Chapman & Hall/CRC, pg 10.

rsync excluding Rdata files

rsync -avz --exclude 'temper.*.Rdata' asthma aitken:~/donnpc/

plot with confidence intervals

## modified version of the code I found on the web.
## I forgot who wrote the origial version

plotCI <- function (x, y = NULL, uiw, liw = uiw, ylim=NULL,..., sfrac = 0.01)  {
 col2 = 'gray'
 if (is.list(x)) {
   y <- x$y
   x <- x$x
 }
 if (is.null(y)) {
   if (is.null(x))
     stop("both x and y NULL")
   y <- as.numeric(x)
   x <- seq(along = x)
 }
 ui <- y + uiw
 li <- y - liw
 if(is.null(ylim)) {
   plot(x, y, ylim = range(c(y, ui, li),na.rm=T), ...)} else
   {plot(x, y, ylim = ylim, ...)}
 smidge <- diff(par("usr")[1:2]) * sfrac
 segments(x, li, x, ui,col=col2)
 x2 <- c(x, x)
 ul <- c(li, ui)
 segments(x2 - smidge, ul, x2 + smidge, ul,col=col2)
 if(is.null(ylim))  ylim=range(c(y, ui, li),na.rm=T)
 points(x, y, ...)
 invisible(list(x = x, y = y))
}

image with gray scale

imagebw = function(x,y,z,...) image(x,y,z,col=gray((33:64)/64),...)

Mac Eye for the Linux Guy

Mac Eye for the Linux Guy

rsync excluding subdirectories contents

rsync -dvlptgoDz source.server:~/source.dir/ target.dir