#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#define RND_MAX 2147480000.0
int main()
{
time_t seconds;
time(&seconds);
srand((unsigned int) seconds);
int i;
float x,y;
for(i=-10;i<=10;i=i+1)
{
/*
if(i%2)
printf("%d 1\n%d 2\n%d 3\n",i,i+1,i+2);
//printf("1\n1\n%d\n",rand()%4);
else
printf("%d 4\n%d 3\n%d 2\n",i,i+1,i+2);
//printf("4\n3\n%d\n",rand()%15);
*/
//printf("%d\n",i*i*i+3*i*i+2);
//printf("%d\n",3*i*i+6*i);
printf("%d\n",6*i+6);
}
return 0;
}
i<- -5:5
x<-i*i*i+3*i*i+5*i-50
fx_spline <- splinefun(x,i)
plot(i,x, type='l' )
#lines(i,fx_spline(x, deriv=1), type='l', col='green')
#lines(i, fx_spline(x, deriv=2), type='l', col='red')
lines(i,diff(x), type='l', col='green')
#lines(i, fx_spline(x, deriv=2), type='l', col='red')
lines(i,3*i*i+6*i+5, type='l', col='blue')
lines(i,6*i+6, type='l', col='orange')
# numerical integration in R
# example based on: http://tolstoy.newcastle.edu.au/R/help/04/10/6138.html
# sample data: FTIR spectra
x <- read.csv(url('http://casoilresource.lawr.ucdavis.edu/drupal/files/fresh_li_material.CSV'), header=FALSE)[100:400,]
names(x) <- c('wavenumber','intensity')
# fit a piece-wise linear function
fx.linear <- approxfun(x$wavenumber, x$intensity)
# integrate this function, over the original limits of x
Fx.linear <- integrate(fx.linear, min(x$wavenumber), max(x$wavenumber))
# fit a smooth spline, and return a function describing it
fx.spline <- splinefun(x$wavenumber, x$intensity)
# integrate this function, over the original limits of x
Fx.spline <- integrate(fx.spline, min(x$wavenumber), max(x$wavenumber))
# visual check, linear and spline fits shifted up for clarity
par(mar=c(0,0,0,0))
plot(x, type = "p", las=1, axes=FALSE, cex=0.5, ylim=c(0,0.12))
lines(x$wavenumber, fx.linear(x$wavenumber) + 0.01, col=2)
lines(x$wavenumber, fx.spline(x$wavenumber) + 0.02, col=3)
grid(nx=10, col=grey(0.5))
legend(x=615, y=0.11, legend=c('original','linear','spline'), col=1:3, pch=c(1,NA,NA), lty=c(NA, 1, 1), bg='white')
# results are pretty close
data.frame(method=c('linear', 'spline'), area=c(Fx.linear$value, Fx.spline$value),error=c(Fx.linear$abs.error,Fx.spline$abs.error))
method area error
1 linear 27.71536 0.0005727738
2 spline 27.71527 0.0030796529
par(mar=c(0,0,0,0), mfcol=c(2,1))
plot(x, type = "l", lwd=2, axes=FALSE)
grid(nx=10, col=grey(0.5))
plot(x$wavenumber, fx.spline(x$wavenumber, deriv=1), type='l', axes=FALSE)
lines(x$wavenumber, fx.spline(x$wavenumber, deriv=2), col='red')
grid(nx=10, col=grey(0.5))
abline(h=0, lty=2)
legend('topright', legend=c('1st derivative','2nd derivative'), lty=1, col=1:2, bg='white')
source : http://casoilresource.lawr.ucdavis.edu/drupal/node/896
features <- function(x, ymat, smoother=c("smooth.spline", "glkerns"), fits.return=FALSE, control=list( ), ...)
# Extracting features of a noisy and discretely sampled function.
# Author: Ravi Varadhan, October, 2010
#
# Function arguments:
# x = vector of independent variable, e.g. time
# ymat = matrix of time-series to be smoothed
# smoother = technique to use for automatic smoothing; only two options are allowed
# fits.return = whether to return the smoothed function and its first and second derivatives
#
# control = a list of control variables describe below
# npts = number of points to use in computing features
# plot.it = a logical variable indicating whether the smoothed results should be plotted
# c.outlier = a constant denoting number of standard deviations away from smooth fit for determining whether a point is an outlier
# decim.out = number of decimals to display in the features output
#
# Extracted features:
# frms = Root-mean-squared function value over the range of data
# fmean = Mean function value over the range of data
# fsd = square-root of the average variance of the function around its mean
# noise = standard deviation of residuals after smoothing
# snr = signal-to-noise ratio = fsd/noise
# fwiggle = root-mean-squared second-derivative of function; a measure of "wiggliness" in the function
# deriv.range = minimum and maximum of first derivative of smoothed function scaled by fsd
# critical.pts = X-locations where the first derivative is zero
# curvature = second derivative of smoothed function at critical points ( > 0 for minimum, and < 0 for maximum)
# scaled by fsd
# outliers = points that are potential outliers
{
ctrl <- list(npts=max(100, 2*length(x)), plot.it=FALSE, c.outlier=3, decim.out=2)
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
npts <- ctrl$npts
plot.it <- ctrl$plot.it
c.outlier <- ctrl$c.outlier
decim.out <- ctrl$decim.out
trapezoid <- function(x,y) sum(diff(x)*(y[-1] + y[-length(y)]))/2
smoother <- match.arg(smoother, c("smooth.spline", "glkerns"))
ipkgs <- installed.packages()
if (smoother == "glkerns") {
if ("lokern" %in% ipkgs[,1]) require(lokern, quietly=TRUE)
else stop("Install package `lokern'", call.=FALSE)
deriv1 <- function(z) glkerns(x, y, deriv=1, x.out=z, hetero=TRUE)$est
deriv2 <- function(z) glkerns(x, y, deriv=2, x.out=z, hetero=TRUE)$est
}
if (smoother == "smooth.spline") {
deriv1 <- function(z) predict(fit, deriv=1, x=z)$y
deriv2 <- function(z) predict(fit, deriv=2, x=z)$y
}
x.out <- seq(min(x), max(x), length=npts)
if (is.null(dim(ymat))) ymat <- matrix(ymat, nrow=1, ncol=length(ymat))
nr <- nrow(ymat)
nc <- ncol(ymat)
fmat <- matrix(NA, nr, 10)
cpmat <- matrix(NA, nr, max(ceiling(nc/10), 100))
cvmat <- matrix(NA, nr, max(ceiling(nc/10), 100))
olmat <- matrix(NA, nr, max(ceiling(nc/10), 100))
if (plot.it & nr > 1) par(ask=TRUE)
if (fits.return) fits <- fits1 <- fits2 <- matrix(NA, nr, nc)
for (k in 1:nr) {
y <- as.numeric(ymat[k, ])
if (smoother == "glkerns") {
fit <- glkerns(x, y, deriv=0, x.out=x, hetero=FALSE)
fit0 <- glkerns(x, y, deriv=0, x.out=x.out, hetero=TRUE)$est
fit1 <- glkerns(x, y, deriv=1, x.out=x.out, hetero=TRUE)$est
fit2 <- glkerns(x, y, deriv=2, x.out=x.out, hetero=TRUE)$est
resid <- y - fit$est
resid.scaled <- abs(scale(resid))
if (fits.return) {
fits[k, ] <- fit$est
fits1[k, ] <- glkerns(x, y, deriv=1, x.out=x, hetero=TRUE)$est
fits2[k, ] <- glkerns(x, y, deriv=2, x.out=x, hetero=TRUE)$est
}
}
if (smoother == "smooth.spline") {
fit <- smooth.spline(x, y, cv=FALSE)
fit0 <- predict(fit, deriv=0, x=x.out)$y
fit1 <- predict(fit, deriv=1, x=x.out)$y
fit2 <- predict(fit, deriv=2, x=x.out)$y
resid <- y - predict(fit)$y
resid.scaled <- abs(scale(resid))
if (fits.return) {
fits[k, ] <- predict(fit)$y
fits1[k, ] <- predict(fit, deriv=1)$y
fits2[k, ] <- predict(fit, deriv=2)$y
}
}
fmean <- trapezoid(x.out, fit0) / diff(range(x.out))
fstar <- sqrt(trapezoid(x.out, fit0^2) / diff(range(x.out)))
fsd <- sqrt(fstar^2 - fmean^2)
noise <- attr(resid.scaled, "scaled:scale")
snr <- fsd/noise
fwiggle <- sqrt(trapezoid(x.out, fit2^2) / diff(range(x.out)))
d0 <- (fit1[-1] * fit1[-npts]) < 0
ncrt <- sum(d0)
if (ncrt == 0) crtpts <- curv <- NULL else {
crtpts <- curv <- rep(NA, ncrt)
ind0 <- (1:npts)[d0]
for (i in 1:ncrt) {
temp <- try(uniroot(interval=c(x.out[ind0[i]], x.out[1 + ind0[i]]), f=deriv1), silent=TRUE)
if (class(temp) != "try-error") {
crtpts[i] <- temp$root
curv[i] <- deriv2(temp$root)
}
}
}
rm(fit)
outl <- resid.scaled > c.outlier
if (sum(outl) == 0) outliers <- NULL else outliers <- x[outl]
#plot
if (plot.it) {
plot(x, y, type="p", ...)
lines(x.out, fit0, col=2)
}
if (!is.null(crtpts) ) {
cpmat[k, 1:length(crtpts)] <- crtpts
cvmat[k, 1:length(crtpts)] <- curv
}
if (!is.null(outliers) ) olmat[k, 1:length(outliers)] <- outliers
fmat[k, ] <- as.numeric(c(fmean, as.numeric(range(fit0)), fsd, noise, snr, as.numeric(range(fit1)), fwiggle, length(crtpts)))
}
fmat <- as.data.frame(fmat)
names(fmat) <- c("fmean", "fmin", "fmax", "fsd", "noise", "snr", "d1min", "d1max", "fwiggle", "ncpts")
ncpmat <- max(apply(cpmat, 1, function(x) sum(!is.na(x))))
ncvmat <- max(apply(cvmat, 1, function(x) sum(!is.na(x))))
nolmat <- max(apply(olmat, 1, function(x) sum(!is.na(x))))
if (ncpmat==0) cpmat <- NULL else cpmat <- round(cpmat[, 1:ncpmat, drop=FALSE], decim.out)
if (ncvmat==0) cvmat <- NULL else cvmat <- round(cvmat[, 1:ncvmat, drop=FALSE], decim.out)
if (nolmat==0) olmat <- NULL else olmat <- olmat[, 1:nolmat, drop=FALSE]
ret.obj <- list(fmat=round(fmat, decim.out), cptmat=cpmat, curvmat=cvmat, olmat=olmat)
if (fits.return) attr(ret.obj, "fits") <- list(x=x, fits0=fits, fits1=fits1, fits2=fits2)
return(ret.obj)
}
---------------------------------------------------------------------------------------------
op <- par(mfrow = c(3, 1))
source("features.R")
#x <- sort(runif(n))
#y <- exp(-0.2 * sin(2*pi*x)) + rnorm(n, sd=0.05)
a<-read.table("WEBM1")
x<- seq(1, 2000, by = 1)
y<-a$V1
#y<-a$V2
#x<- seq(-5, 5, by = 0.01)
#y<- x*x*x
ans1 <- features(x, y, fits.return=TRUE, control=list(plot.it=FALSE))
#ans <- features(x, y, fits.return=TRUE, control=list(plot.it=TRUE))
fits1 <- attr(ans1, "fits")
#yexact <- -0.2 * 2*pi * cos(2*pi*x) * exp(-0.2 * sin(2*pi*x))
#lines(x, yexact, col=3)
plot(x,y,type="l",col=2,ylab="Original data",main="Video Trace")
grid(col=grey(0.3))
plot(x, fits1$fits1, type="l", xlab="x",col=3, ylab="First derivative")
#plot(fits1$x, fits1$fits1, type="l", xlab="x",col=3, ylab="First derivative")
grid(col=grey(0.3))
ans2 <- features(x, fits1$fits1, fits.return=TRUE, control=list(plot.it=FALSE))
fits2 <- attr(ans2, "fits")
plot(x, fits2$fits1, type="l", xlab="x",col=3, ylab="2nd derivative")
grid(col=grey(0.3))
ans3 <- features(x, fits2$fits1, fits.return=TRUE, control=list(plot.it=FALSE))
fits3 <- attr(ans3, "fits")
#plot(x, fits3$fits1, type="l", xlab="x",col=3, ylab="3rd derivative")
#grid(col=grey(0.3))
par(op)
-------------------------------------------------------------------------------------
op <- par(mfrow = c(2, 1))
x<-scan("a15_m4_mpeg")
y<- seq(1, 299, by = 1)
source("features.R")
ans1 <- features(x, y, fits.return=TRUE, control=list(plot.it=FALSE))
plot(x,y,type="l",col=2,ylab="Original data")
grid(col=grey(0.3))
par(op)
-------------------------------------------------------------------------------------
op <- par(mfrow = c(3, 1))
x <- seq(0, 9, by = 0.1)
# u is the true function, u1-u3 are its symbolic derivatives wrt x
# u1 = f', u2 = f'', u3 = f'''
#u <- expression(x*x*x+3*x*x+5*x-50)
u <- expression(sin(pi * (x - 0.5)))
u1 <- as.expression(D(u, 'x'))
u2 <- as.expression(D(u1, 'x'))
u3 <- as.expression(D(u2, 'x'))
#plot( eval(u), type = 'l',main = expression(f(x) == sin(pi * (x - 0.5))))
plot(eval(u), type = 'l', ylim = c(-30, 30),ylab = "", xlab = "",main = expression(f(x) == sin(pi * (x - 0.5))))
lines( eval(u1), type = 'l', col = 'red')
lines(eval(u2), type = 'l', col = 'blue')
lines( eval(u3), type = 'l', col = 'green')
grid(col=grey(0.3))
#legend('topright', legend = c("f", "f'", "f''", "f'''"),col = c('black', 'red', 'blue', 'green'), lwd = 2)
# y2 is an evaluation of u at each x point
y2 <- sin(pi * (x - 0.5))
#y2 <- x*x*x+3*x*x+5*x-50
# Construct the spline function approximation to u
f <- splinefun(x, y2)
# Plot the interpolation function and its first three derivatives
curve(f(x), 1, 10, col = "black", lwd = 1.5, ylim = c(-30, 30), xlim = c(2,11),ylab = "", main = 'Spline interpolation of f')
curve(f(x, deriv=1), 1, 10, col = 'red', lwd = 1.5, add = TRUE)
curve(f(x, deriv=2), 1, 10, col = 'blue', lwd = 1.5, add = TRUE)
curve(f(x, deriv=3), 1, 10, col = 'green', lwd = 1.5, add = TRUE)
grid(col=grey(0.3))
x<- seq(0, 9, by = 0.01)
plot(sin(pi * (x - 0.5)), type = 'l', ylim = c(-30, 30),ylab = "", xlab = "",main ='Direct Derivatives')
lines((cos(pi * (x - 0.5))*pi), type='l', col='red')
lines((-1)*(sin(pi * (x - 0.5))*(pi*pi)), type='l', col='blue')
lines((-1)*(cos(pi * (x - 0.5))*(pi*pi*pi)), type='l', col='green')
grid(col=grey(0.3))
par(op)