Hi,
If you didn't receive the attachment properly, here it is again.
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Ravi Varadhan
Sent: Wednesday, March 05, 2008 4:58 PM
To: 'Spencer Graves'; 'Levi Waldron'
Cc: 'R-help mailing list'
Subject: Re: [R] differentiating a numeric vector
Hi,
Here is another approach, in addition to the suggestions made by Spencer and
Gabor. It uses the spm() function in SemiPar package. An advantage of this
approach is that the smoothing parameter is automatically estimated using
REML (here I use default knot locations, but this can be specified
explicitly).
You also need to source in the plot.spm() function that I have created by
slightly modifying the plot.spm() that comes with SemiPar. This is
necessary to obtain numerical values of derivatives at x-locations in
addition to simply plotting the derivative curves.
library(SemiPar)
source("plotspm.r")
# An example
k <- 10
x <- sqrt(runif(500))
y <- pnorm(x) + sin(k*pi*x^2) + rnorm(500,mean=0,sd=0.5)
fit<-spm(y ~ f(x, basis="trunc.poly", degree=3), omit.missing=TRUE)
deriv <- plot(fit, drv=1) # plot and store first derivative
deriv2 <- plot(fit, drv=2) # plot and store second derivative
Hope this helps,
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Spencer Graves
Sent: Wednesday, March 05, 2008 4:11 PM
To: Levi Waldron
Cc: R-help mailing list
Subject: Re: [R] differentiating a numeric vector
Have you looked at the 'fda' package? It has many functions for
doing what you want. A strength is that it is a companion package for
two books on that and related issues, and includes script files under
"~R.installation.directory\library\fda\scripts" to reproduce some of the
analyses. It may include more than you want to consider, but for me,
too much is usually better than nothing.
hope this helps.
Spencer Graves
p.s. If you try it and have trouble, please submit another question
including commented, minimal, self-contained, reproducible code, as
requested in the posting guide
http://www.R-project.org/posting-guide.html.
Levi Waldron wrote:
> What functions exist for differentiating a numeric vector (in my case
> spectral data)? That is, experimental data without an analytical
> function. ie,
>
>
>> x <- seq(1,10,0.1)
>> y=x^3+rnorm(length(x),sd=0.01) #although the real function would be
nothing simple like x^3...
>> derivy <- ....
>>
>
> I know I could just use diff(y) but it would be nice to estimate
> derivatives at the endpoints, calculate higher-order derivatives, and
> maybe have some smoothing options ie by differentiating a spline or
> something like that.
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
plot.spm <- function (x, ...)
{
object <- x
plot.params <- plotControl(...)
drv <- plot.params$drv
se <- plot.params$se
shade <- plot.params$shade
default.ylim <- FALSE
if (is.null(plot.params$ylim))
default.ylim <- TRUE
lin.present <- !is.null(object$info$lin)
pen.present <- !is.null(object$info$pen)
krige.present <- !is.null(object$info$krige)
random.present <- !is.null(object$info$random)
const.only <- ((!lin.present) & (!pen.present) & (!krige.present) &
(!random.present))
block.inds <- object$aux$block.inds
if (pen.present | krige.present)
coefs <- c(object$fit$coef$fixed, object$fit$coef$random)
else coefs <- object$fit$coef$fixed
if (se == TRUE)
cov.mat <- object$aux$cov.mat
if (random.present) {
random.inds <- block.inds[[length(block.inds)]]
coefs <- coefs[-random.inds]
if (se == TRUE)
cov.mat <- cov.mat[-random.inds, -random.inds]
}
basis <- object$info$pen$basis
if (drv == 0)
ave.vals <- 1
if (drv > 0)
ave.vals <- 0
if (lin.present)
ave.vals <- c(ave.vals, apply(object$info$lin$x, 2, mean))
if (pen.present) {
num.pen <- ncol(as.matrix(object$info$pen$x))
for (ipen in 1:num.pen) {
deg.val <- object$info$pen$degree[ipen]
if (basis == "trunc.poly")
ncol.X.val <- deg.val
if (basis == "tps")
ncol.X.val <- (deg.val - 1)/2
ave.val <- mean(as.matrix(object$info$pen$x)[, ipen])
ave.vals <- c(ave.vals, ave.val)
if (ncol.X.val > 1)
for (ipow in 2:ncol.X.val) ave.vals <- c(ave.vals,
ave.val^ipow)
}
}
if (!pen.present)
num.pen <- 0
if (krige.present) {
x1.v <- object$info$krige$x[, 1]
x2.v <- object$info$krige$x[, 2]
m.krige <- (object$info$krige$degree/2) + 1
X.kg <- NULL
for (im in 2:m.krige) {
X.kg <- cbind(X.kg, x1.v^(im - 1))
if (im > 2)
for (ipow in 2:(im - 1)) X.kg <- cbind(X.kg,
(x1.v^(im - ipow)) * (x2.v^(ipow - 1)))
X.kg <- cbind(X.kg, x2.v^(im - 1))
}
ave.vals <- c(ave.vals, apply(X.kg, 2, mean))
}
if (pen.present) {
num.pen <- ncol(as.matrix(object$info$pen$x))
for (ipen in 1:num.pen) {
deg.val <- object$info$pen$degree[ipen]
knots <- object$info$pen$knots[[ipen]]
ave.val <- mean(as.matrix(object$info$pen$x)[, ipen])
if (basis == "trunc.poly") {
ave.val.knots <- ave.val - knots
ave.val.knots <- (ave.val.knots * (ave.val.knots >
0))^deg.val
}
if (basis == "tps") {
ave.val.knots <- abs(ave.val - knots)^deg.val
sqrt.Omega <- object$info$trans.mat[[ipen]]
ave.val.knots <- t(solve(sqrt.Omega, ave.val.knots))
}
ave.vals <- c(ave.vals, ave.val.knots)
}
}
if (krige.present) {
ave.val.x1 <- mean(object$info$krige$x[, 1])
ave.val.x2 <- mean(object$info$krige$x[, 2])
knots <- object$info$krige$knots
knots.1 <- knots[, 1]
knots.2 <- knots[, 2]
num.knots <- length(knots.1)
ave.val.knots <- NULL
dists.ave <- sqrt((ave.val.x1 - knots.1)^2 + (ave.val.x2 -
knots.2)^2)
ave.val.knots <- tps.cov(dists.ave, m = m.krige, d = 2)
ave.val.knots <-
solve(object$info$trans.mat[[length(object$info$trans.mat)]],
as.matrix(ave.val.knots))
ave.vals <- c(ave.vals, ave.val.knots)
}
if (const.only) {
mean.est <- coefs[1]
x.grid <- c(0.25, 0.75)
se.val <- sqrt(diag(object$aux$cov.mat))[1]
lower <- mean.est - 2 * se.val
upper <- mean.est + 2 * se.val
plot(0, 0, type = "n", xlim = c(0, 1), ylim = c(mean.est -
2.5 * se.val, mean.est + 2.5 * se.val), xaxt = "n",
xlab = "", ylab = "", bty = "l")
if (se == TRUE) {
if (shade == FALSE) {
lines(x.grid, rep(lower, 2), lty = 2, lwd = 2,
col = "black")
lines(x.grid, rep(upper, 2), lty = 2, lwd = 2,
col = "black")
}
if (shade == TRUE)
polygon(c(x.grid, rev(x.grid)), c(rep(lower,
2), rep(upper, 2)), border = FALSE, col = "grey70")
}
lines(x.grid, rep(mean.est, 2), lwd = 2)
}
if (lin.present | pen.present) {
curve.type <- NULL
if (lin.present)
curve.type <- c(curve.type, rep("lin",
ncol(as.matrix(object$info$lin$x))))
if (pen.present)
curve.type <- c(curve.type, rep("pen",
ncol(as.matrix(object$info$pen$x))))
num.curves <- length(curve.type)
}
if (!(lin.present | pen.present))
num.curves <- 0
plot.params <- set.plot.dflts(object, plot.params, num.curves)
if ((!is.null(plot.params$xlim)) & (!is.list(plot.params$xlim))) {
if (length(plot.params$xlim) != 2)
stop("illegal xlim")
plot.params$xlim <- list(lower = rep(plot.params$xlim[1],
num.curves), upper = rep(plot.params$xlim[2], num.curves))
}
if ((!is.null(plot.params$ylim)) & (!default.ylim) &
(!is.list(plot.params$ylim))) {
if (length(plot.params$ylim) != 2)
stop("illegal ylim")
plot.params$ylim <- list(lower = rep(plot.params$ylim[1],
num.curves), upper = rep(plot.params$ylim[2], num.curves))
}
if (lin.present | pen.present) {
x.vals <- NULL
if (lin.present)
x.vals <- cbind(x.vals, object$info$lin$x)
if (pen.present)
x.vals <- cbind(x.vals, object$info$pen$x)
fc.stt.pos <- 2
pen.num <- 1
for (j in 1:num.curves) {
grid.size <- plot.params$grid.size[j]
if (drv == 0)
C.grid <- matrix(rep(ave.vals, grid.size), grid.size,
length(ave.vals), byrow = TRUE)
if (drv > 0)
C.grid <- matrix(0, grid.size, length(ave.vals))
x.grid <- seq(plot.params$xlim$lower[j], plot.params$xlim$upper[j],
length = plot.params$grid.size[j])
if (curve.type[j] == "lin")
deg.val <- 1
if (curve.type[j] == "pen")
deg.val <- object$info$pen$degree[pen.num]
X.g.inst <- NULL
if (curve.type[j] == "lin") {
if (drv == 0)
X.g.inst <- cbind(X.g.inst, x.grid)
if (drv == 1)
X.g.inst <- cbind(X.g.inst, rep(1, grid.size[j]))
if (drv > 1)
X.g.inst <- cbind(X.g.inst, rep(0, grid.size[j]))
}
if (curve.type[j] == "pen") {
if (basis == "trunc.poly")
ncol.X.val <- deg.val
if (basis == "tps")
ncol.X.val <- (deg.val - 1)/2
if (drv == 0) {
for (pow in 1:ncol.X.val) {
new.col.data <- x.vals[, j]^pow
new.col <- x.grid^pow
X.g.inst <- cbind(X.g.inst, new.col)
}
}
if (drv > 0) {
for (pow in 1:ncol.X.val) {
new.col.data <- x.vals[, j]^pow
pow.drv <- pow - drv
if (pow.drv >= 0)
new.col <- prod(pow:(pow.drv + 1)) * (x.grid^pow.drv)
else new.col <- rep(0, length(x.grid))
X.g.inst <- cbind(X.g.inst, new.col)
}
}
}
C.g.inst <- X.g.inst
if (curve.type[j] == "lin") {
inst.col.inds <- fc.stt.pos
fc.stt.pos <- fc.stt.pos + 1
}
if (curve.type[j] == "pen") {
if (basis == "trunc.poly")
ncol.X.val <- deg.val
if (basis == "tps")
ncol.X.val <- (deg.val - 1)/2
inst.col.inds <- fc.stt.pos:(fc.stt.pos + ncol.X.val -
1)
fc.stt.pos <- fc.stt.pos + ncol.X.val
knots <- object$info$pen$knots[[pen.num]]
num.knots <- length(knots)
Z.g.inst <- outer(x.grid, knots, "-")
if (basis == "trunc.poly") {
if (drv == 0)
Z.g.inst <- (Z.g.inst * (Z.g.inst > 0))^deg.val
else {
if (deg.val >= drv) {
mfac <- prod((deg.val - drv + 1):deg.val)
Z.g.inst <- mfac * (Z.g.inst * (Z.g.inst >
0))^(deg.val - drv)
}
else Z.g.inst <- matrix(0, length(x.grid),
length(knots))
}
}
if (basis == "tps") {
if (drv == 0) {
Z.g.inst <- abs(Z.g.inst)^deg.val
sqrt.Omega <- object$info$trans.mat[[pen.num]]
Z.g.inst <- t(solve(sqrt.Omega, t(Z.g.inst)))
}
else {
if (deg.val >= drv) {
mfac <- prod((deg.val - drv + 1):deg.val)
Z.g.inst <- mfac * (Z.g.inst^(deg.val -
drv - 1) * abs(Z.g.inst))
sqrt.Omega <- object$info$trans.mat[[pen.num]]
Z.g.inst <- t(solve(sqrt.Omega, t(Z.g.inst)))
}
else Z.g.inst <- matrix(0, length(x.grid),
length(knots))
}
}
C.g.inst <- cbind(C.g.inst, Z.g.inst)
inst.col.inds <- c(inst.col.inds, block.inds[[j +
1]][-(1:ncol.X.val)])
pen.num <- pen.num + 1
}
C.grid[, inst.col.inds] <- C.g.inst
y.grid <- as.vector(C.grid %*% coefs)
if (default.ylim) {
plot.params$ylim$lower[j] <- min(y.grid)
plot.params$ylim$upper[j] <- max(y.grid)
}
if (se == TRUE) {
se.grid <- sqrt(diag(C.grid %*% cov.mat %*% t(C.grid)))
lower <- y.grid - 2 * se.grid
upper <- y.grid + 2 * se.grid
if (default.ylim) {
plot.params$ylim$lower[j] <- min(lower)
plot.params$ylim$upper[j] <- max(upper)
}
ylim.val <- range(c(lower, upper))
}
if (plot.params$plot.it[j] == TRUE) {
plot(x.grid, y.grid, type = "n", bty = plot.params$bty[j],
main = plot.params$main[j], xlab = plot.params$xlab[j],
ylab = plot.params$ylab[j], xlim =
c(plot.params$xlim$lower[j],
plot.params$xlim$upper[j]), ylim =
c(plot.params$ylim$lower[j],
plot.params$ylim$upper[j]))
if (se == TRUE) {
if (shade == FALSE) {
lines(x.grid, lower, lty = plot.params$se.lty[j],
lwd = plot.params$se.lwd[j], col = plot.params$se.col[j])
lines(x.grid, upper, lty = plot.params$se.lty[j],
lwd = plot.params$se.lwd[j], col = plot.params$se.col[j])
}
if (shade == TRUE) {
mat <- cbind(x.grid, lower, upper)
mat <- mat[sort.list(mat[, 1]), ]
x.grid <- mat[, 1]
lower <- mat[, 2]
upper <- mat[, 3]
polygon(c(x.grid, rev(x.grid)), c(lower,
rev(upper)), col = plot.params$shade.col[j],
border = FALSE)
}
}
if ((drv >= 1) & (plot.params$zero.line == TRUE))
abline(h = 0, err = -1)
lines(x.grid, y.grid, lty = plot.params$lty[j],
lwd = plot.params$lwd[j], col = plot.params$col[j])
if (plot.params$jitter.rug == TRUE)
rug.x <- jitter(x.vals[, j])
if (plot.params$jitter.rug == FALSE)
rug.x <- x.vals[, j]
rug(rug.x, quiet = 1, col = plot.params$rug.col[j])
}
}
}
if (krige.present) {
if (plot.params$plot.image == TRUE) {
if ((!lin.present) & (!pen.present))
num.curves <- 0
pixel.info <- get.pixel.info(plot.params)
on.inds <- pixel.info$on.inds
off.inds <- pixel.info$off.inds
image.grid.size <- plot.params$image.grid.size
num.pix <- image.grid.size[1] * image.grid.size[2]
num.on.pix <- length(on.inds)
C.grid <- matrix(rep(ave.vals, num.on.pix), num.on.pix,
length(ave.vals), byrow = TRUE)
z <- rep(0, num.pix)
X.g.inst <- as.matrix(pixel.info$newdata)
x1.g.inst <- X.g.inst[, 1]
x2.g.inst <- X.g.inst[, 2]
m.krige <- (object$info$krige$degree/2) + 1
if (m.krige > 2)
for (im in 3:m.krige) {
X.g.inst <- cbind(X.g.inst, x1.g.inst^(im -
1))
if (im > 2)
for (ipow in 2:(im - 1)) X.g.inst <- cbind(X.g.inst,
(x1.g.inst^(im - ipow)) * (x2.g.inst^(ipow -
1)))
X.g.inst <- cbind(X.g.inst, x2.g.inst^(im -
1))
}
knots <- object$info$krige$knots
knots.1 <- knots[, 1]
knots.2 <- knots[, 2]
num.knots <- length(knots.1)
dists.g <- outer(x1.g.inst, knots.1, "-")^2
dists.g <- sqrt(dists.g + outer(x2.g.inst, knots.2,
"-")^2)
prelim.Z.g.inst <- tps.cov(dists.g, m = m.krige,
d = 2)
Z.g.inst <- t(solve(object$info$trans.mat[[num.pen +
1]], t(prelim.Z.g.inst)))
C.g.inst <- cbind(X.g.inst, Z.g.inst)
inst.col.inds <- block.inds[[num.curves + 2]]
C.grid[, inst.col.inds] <- C.g.inst
z[on.inds] <- as.vector(C.grid %*% coefs)
z[off.inds] <- NA
z <- matrix(z, image.grid.size[1], image.grid.size[2])
imageFull(z, plot.params)
}
}
return(list(xval=x.grid,yval=y.grid,se=se.grid))
invisible()
}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
