**This is for testing various techniques:**

Upload of an equation from a pdf converted to a jpg

Seems to work reasonably well. Now a slightly smaller png file:

Not bad.

Try WordPress latex:

Now redo the equation above:

Yes!! It works.

** **

Post an R script as a block quote:

download.file(“http://www.climateaudit.org/data/antarctic_mann/Data.tab”,”temp2.dat”,mode=”wb”);load(“temp2.dat”)

slopes = function(tsdat) {

star = min(time(tsdat))

slopes = rep(NA,12)

for (i in 1:12){

dat.win = window(tsdat,start=c(star,i),deltat=1)

yr = time(dat.win)

slopes[i] = ifelse(sum(!is.na(dat.win)) < 2, NA, coef(lm(dat.win~yr))[2])

}

slopes}

slope(Data$aws[,8])

all.slope = function(alldat) {

nc = ncol(alldat)

slop = matrix(NA, nrow = nc, ncol = 12)

colnames(slop) = month.abb

for (j in 1:nc) {

slop[j,]= slope(alldat[,j])}

slop}

test3 = all.slope(Data$aws[,6])

Quotes copy out wrong. This needs fixing.

New test:

Attach text document as .doc file:

collation_functions

### Like this:

Like Loading...

*Related*

Roman,

I have written a R script to verify Giorgio’s calculation. It does – I get a symmetric distribution, mean 0.0175 deg C/decade. and standard deviation 0.189 C/dec. The histogram is here. I tried to post a comment on his blog, but it hasn’t shown yet.

To make it easier in R I edited (with emacs) to turn -9999 into NA and to separate the station numbers from the years.

Here’s the script:

#### A program written by Nick Stokes, 13 Dec 2009, to calculate the changes to regression

# slopes caused by adjustments to the GHCN temperatures v2.mean_adj-v2.mean

# A function to calculate regression slope. I hope it is faster than lm()

slope<-function(v,jj){

m=jj-mean(jj)

s=(v %*% m)/(m %*% m)

s

}

#####################

# read data from v2.mean and v2.mean_adj, downloaded from http://www1.ncdc.noaa.gov/pub/data/ghcn/v2/

# I edited (emacs) to put a blank between the station number and year, and to change -9999 to NA (add .txt)

# Read in data from the files in matrix form

if(T){ #change to F after you have read in th efiles once

vmean <- matrix(scan("v2.mean.txt", 0, skip=0,na.strings = "NA"), ncol=14, byrow=TRUE)

vmean_adj <- matrix(scan("v2.mean_adj.txt", 0, skip=0,na.strings = "NA"), ncol=14, byrow=TRUE)

# Now, to save time, move to annual averages

vmean_ann=vmean[,1:3]

vmean_ann[,3]=rowMeans(vmean[,3:14], na.rm = T)

vmean_ann_adj=vmean_adj[,1:3]

vmean_ann_adj[,3]=rowMeans(vmean_adj[,3:14], na.rm = T)

}

# Initialise

vv=rep(0.,200) # regression y vector

jj=rep(0,200) # regression y vector

grad=rep(0.,9999) # gradients (the output result)

len=length(vmean_ann[,1])

jmax=length(vmean_ann_adj[,1])

j=1

k=0

kk=0

m=0

# counters. j is row of adjusted file. m is station counter

# k,kk are local row (year) counter (for station m). k skips NA's, kk doesn't

# loop over all rows in v2.mean

for(i in 1:(len-1)){

kk=kk+1

# to find matching rows, first check diff between stat nos and years

u=vmean_ann_adj[j,]-vmean_ann[i,]

# If the adjusted counter has got ahead of the unadj, wait

if(u[1]<0){

if(j<jmax)j=j+1; u=vmean_ann_adj[j,]-vmean_ann[i,]

}

# If we have a match, add to regression vec vv[]

if(u[1]==0 & u[2]==0 ){

if(!is.na(u[3])){ # don't add to regression if NA

k=k+1 # local adjusted counter

jj[k]=kk # x for regression

vv[k]=u[3] # discrepancies for regression

}

if(j0){

m=m+1 # m is station counter

grad[m]=slope(vv[1:k],jj[1:k]) # compute regression slope

k=0 # zero local counters

kk=0

}

}

# Now prepare histogram. Comment out jpeg and dev.off() to get screen graphics

jpeg(“GHCNAdjustments.jpg”)

hist(grad[1:m],nclass=200,xlab=”degrees C/decade”,main=”GHCN adjustment change to trend”) # draw histogram

a=c(mean(grad[1:m])); a # Mean slope change

dev.off()

Nick, I haven’t been ignoring your fine R script work. It is just that this post is in a side thread which I have been using to test features of HTML, etc. and was not intended to be used for comments.

I haven’t figured out how to move it to the relevant thread yet.