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

Thursday, December 21, 2006

R documentation with graphics



R documentation with graphics

specifying printer lpr


lpr -#2 -sP dj thesis.txt
This command will create a symbolic link to the file thesis.txt in the spool
directory for the printer named dj, where it would be processed by lpd.
It would then print a second copy of thesis.txt.

Monday, December 18, 2006

bold greek letters latex



% use bm package
\usepackage{bm}

% inside math mode
\[ ... {\bm beta} ... \]

better than \boldmath because it works within math mode directly,
no need to use in text mode.

Friday, December 15, 2006

logit - expit


%% maps (0,b) into (-inf,inf)
%% inverse of expit
%% usage: logit(x,b) or logit(x)
%% b defaults to 1
function y = logit(x,b)
if nargin<2
b = 1;
end
y = log( x ./ (b - x) );

%% maps (-inf, inf) into (0,b)
%% inverse of logit
%% usage: expit(x,b) or expit(x)
%% b defaults to 1
function y = expit(x,b)
if nargin<2
b = 1;
end
y = b./(1+exp(-x)) ;

Monday, December 11, 2006

xcdroast

to copy files to CD in linux
> xcdroast

Tuesday, December 5, 2006

psnup

psnup -d -8 oldfile.ps newfile.ps
many postscript pages in one

Friday, November 17, 2006

R equivalent of repmat (matlab)

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

Wednesday, December 8, 2004

ssh without password howto

 
1. Run "ssh-keygen -t rsa".

2. Save the generated key in the default location.

3. Leave the passphrase blank.

4. Change to the .ssh directory, and append the contents of id_rsa.pub
to the file "authorized_keys." If no authorized_keys file exists, you
can just run "cp id_rsa.pub authorized_keys".

Now you should be able to ssh to department machines without typing
your password. Note that if you're doing this from a non-department
Linux machine, you should do steps 1, 2, and 3 from the command prompt
of your machine, and step 4 from the command prompt of any department
Linux machine. If you're doing this from a department Linux machine,
do all steps on that machine.

- Elliot

R legacy help galton

## calling C functions from R
gcc -c kendall.c
R CMD SHLIB -o libkendall.so kendall.o
R
dyn.load('libkendall.so')
xd <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1)
yd <- c( 2.6,  3.1,  2.5,  5.0,  3.6,  4.0,  5.2,  2.8,  3.8)
.C("kendall",as.double(xd),as.double(yd),as.integer(9),
result=as.integer(1))$result

## header
rm(list=ls())
options(contrasts=c("contr.treatment","contr.poly"))
setwd("/mnt/win_c/my documents/stat/ammonia")

## plot function
plot(fun,-5, 5, ylim = c(-.2, 1))

eval.parent(expr) 
## example eval.parent(media <- 1), used inside
## functions to change global variables
media <<- 1 ## does the same

rep(x,n) replicates x n-times

mtext("texto",side=3,line=.5)

install.packages("fields","/aa/haky/R/Library")
library("fields",lib="/aa/haky/R/Library")
detach(package:fields)

setenv R_LIBS /aa/haky/R/Library (then no need to specify lib=...)

xmat <- matrix(NA,n,n)

#print
dev.print()
dev.copy2eps(file="filename.ps")
#
postscript(file="plot.ps",horizontal=F,heigth=9,width=7)
dev.off()

as.numeric(x<1)

# stop watch
tic<-Sys.time()
....do something here
toc<-Sys.time()
print(toc-tic)

paste(s1,s2,sep="")

system("mail -s something haky@galton.uchicago.edu < output")<

## strip an object's attributes:
attributes(x) <- NULL
x # now just a vector of length 6

kk[names(kk)=="site"] is same as kk$site

unlist(x)

debug(myfunction)
undebug(myfunction)

call R from emacs (ESS): alt-x R

par(mfrow=c(3,2))

legacy galton help

# synchronize file system
rsync -avz --exclude '*.eps' Penalty aitken:~/donnpc/ .
# dry run - shows only what would have been transferred
rsync -avzn --exclude '*.eps' Penalty aitken:~/donnpc/ . | less

# indent region mode dependent
M - C - \

# rigid tab region
C - x - TAB

# restart pcmcia, network
/etc/rc.d/init.d/pcmcia restart
/etc/init.d/network restart

# create symbolic links
ln -s /home/azubrow/CMAQ_v4.3/data/emis/M_36_EUSA ./emisfiles

# run identical command in all nodes
sshall uptime
sshall w

To run batch jobs in ravana
chmod u+x batch
./batch &

batch:
#!/bin/sh
/usr/bin/nohup nice -5 R BATCH mplefft.r output

Ctrl-z (suspend)
bg (sends process to background)
fg(bring process to foreground)

rm
rmdir
mkdir
mv
chmod 700 mydir (no access to others)

ftp
$ scp * ravana:~/

gunzip filename: unzip
gzip filename: zip
tar xvf filename.tar: extract tar
tar cvf filename.tar filename or directory: merges files

get rid of the ^M:  in vi type
:%s/(ctrl-v)(ctrl-m)//g
in emacs
M-x replace-regexp RET C-q C-m $ RET RET

chsh #changes shell in unix

mail -s "subject" email@galton <> mkdir work_dir
> cd work_dir
> Splus5 CHAPTER
> touch .Data/.Audit
> chmod 0 .Data/.Audit
emacs kk.s   alt-x S+6

Wednesday, June 2, 2004

cdms quick tutorial

> cdat This starts the CDAT shell, identical to the Python shell.
> import cdms Imports the CDMS module.
> f=cdms.open('wind_comps.nc') Opens the NetCDF file and assigns the file object to variable 'f'.
> f.listvariables() Shows the variable names in the file.
['u', 'v']
> u_wind=f('u') Reads the variable 'u' and assigns it to the Python variable 'u_wind'.
> u_wind.attributes Displays all the attributes pertaining to the variable 'u_wind'.
> v_wind=f('v') Reads the variable 'v' and assigns it to the Python variable 'v_wind'.
> wspd=(u_wind**2+v_wind**2)**0.5
Calculates the combined wind speed and assigns it to the Python variable 'wspd'.
> wspd.id='wspd' Sets the 'id' attribute on the 'wspd' variable.
> wspd.long_name='Wind speed' Sets the 'long_name' attribute on the 'wspd' variable.
> wspd.units='m s**-1' Sets the 'units' attribute on the 'wspd' variable.
> import vcs Imports the VCS (Visualisation Control System) module.
> p=vcs.init() Initialise a VCS canvas and assign Python variable 'p' to it.
> p.plot(wspd) Plot the wind speed variable on canvas 'p'.
> fout=cdms.open('output.nc', 'w') Open a file to write out the wind speed variable to.
> fout.write(wspd) Write the wind speed variable to the output file.

> fout.close() Close the output file.
> CTRL^D Press Control and D to close CDAT.

Contributors

google