Fitting a (non-linear) curve to data points

Download data files here.
Download R script here.

library(chron)
url <- paste("http://www1.ncdc.noaa.gov/pub/orders/cdo/294541.dat", sep="")

dat <- read.table(url, sep = "", skip = 2, header = F)
bostonTemp <- data.frame(dat[c("V11", "V13")])
headers <- read.table(url, colClasses = "character", sep = "",
                      nrow = 1, header = F)
names(bostonTemp) <- c(headers[c("V6", "V8")])
bostonTemp$TMAX <- bostonTemp$TMAX/10
bostonTemp$DATE <- as.chron(paste(substr(bostonTemp$DATE, 5,6),
                                  substr(bostonTemp$DATE, 7,8),
                                  substr(bostonTemp$DATE, 1,4), sep="/"))

par(font=2, font.axis = 2, tck = 0.01)
plot(bostonTemp$DATE, bostonTemp$TMAX)

bostonTemp2013 <- subset(bostonTemp, 
                         bostonTemp$DATE >= as.chron("01/01/2013") & 
                           bostonTemp$DATE <= as.chron("12/31/2013"))
plot(bostonTemp2013$DATE, bostonTemp2013$TMAX) 

x <- as.numeric(bostonTemp2013$DATE)
y <- bostonTemp2013$TMAX
h1 <- mean(x)
k1 <- max(y)
a1 <- -1

yfit<- nls(y ~ a*(x-h)^2 + k, start = list(a = a1, h = h1, k = k1))
yfitParm <- summary(yfit)$para[,1]

ymod <- yfitParm[1]*(x-yfitParm[2])^2 + yfitParm[3]

par(font=2, font.axis = 2, tck = 0.01)
plot(bostonTemp2013$DATE, bostonTemp2013$TMAX)
lines(x, ymod)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s