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

1 comment:

irtaza said...

This article is very helpful many sites provide free online courses with certificates of completion.

Contributors

google