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)