Friday, July 13, 2007

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)
}

Wednesday, July 11, 2007

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)

Monday, June 25, 2007

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 ) )

Monday, June 18, 2007

theoretical ACF of ARMA

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

Friday, June 8, 2007

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 

Tuesday, May 15, 2007

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)}

Thursday, May 3, 2007

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.

Wednesday, March 28, 2007

rsync excluding Rdata files


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

Friday, March 2, 2007

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))
}

Wednesday, February 28, 2007

image with gray scale

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

Thursday, February 15, 2007

Contributors

google