[Swift-commit] r2853 - in SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2: . peak_scripts
noreply at svn.ci.uchicago.edu
noreply at svn.ci.uchicago.edu
Thu Apr 9 16:37:19 CDT 2009
Author: andric
Date: 2009-04-09 16:37:19 -0500 (Thu, 09 Apr 2009)
New Revision: 2853
Added:
SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/
SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/Shellpeak.R
SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/Shellpeak_PK2runbyrun.R
SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/peakfit_2.1.R
SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/preprocess.Enewsmooth.R
Log:
scripts dir for doing peak analysis
Added: SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/Shellpeak.R
===================================================================
--- SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/Shellpeak.R (rev 0)
+++ SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/Shellpeak.R 2009-04-09 21:37:19 UTC (rev 2853)
@@ -0,0 +1,59 @@
+library(stats)
+source("scripts/preprocess.Enewsmooth.R")
+source("scripts/peakfit_2.1.R")
+
+# get quick and dirty estimate for where baseline should be for the assessment.
+# because the cycle is 17on,10off, 63% of points should be above whatever baseline really is
+# meaning that if we sort the y-=values of the 108TR series, it should be rank #33
+
+allinputs <- Sys.getenv("R_SWIFT_ARGS")
+inputfilename <- noquote(strsplit(allinputs," ")[[1]][1])
+outputPDFfilename <- noquote(strsplit(allinputs," ")[[1]][2])
+outputAfilename <- noquote(strsplit(allinputs," ")[[1]][3])
+outputBfilename <- noquote(strsplit(allinputs," ")[[1]][4])
+
+
+dd <- read.table(inputfilename)
+## make sure to change line below to fit the number of timepoints in ts
+quickbase <- (sort(dd[,1])[66])
+adjusted.input = dd[,1]-quickbase
+res_preprocess <- preprocess_1(input = adjusted.input);
+
+test2_res_D2_nlminb <- nlminb(
+start = res_preprocess$guess,
+objective = fcn_L2,
+gradient = test_D_fcn_L2,
+#hessian = test_D2_fcn_L2,
+control = list(eval.max = 50000, iter.max = 5000),
+# the lower line with fewer iterations is just for testing purposes. the upper line is the real one
+# control = list(eval.max = 100, iter.max = 50),
+lower = rep(c(0),length(res_preprocess$guess)),
+fcall = f,
+n = res_preprocess$n,
+x = res_preprocess$interpolated_x,
+y_and_noise = res_preprocess$interpolated_y)
+
+pdf(file=outputPDFfilename)
+op <- par(mfrow=c(1,1))
+plot(res_preprocess$interpolated_x,res_preprocess$interpolated_y,xlab = "",ylab = "BOLD (1E-3)",type = "b");
+plot.est(par = test2_res_D2_nlminb$par, fcall = f, n = res_preprocess$n, res_preprocess$interpolated_x,col = "blue");
+plot.est.individual(par = test2_res_D2_nlminb$par, fcall = f, n = res_preprocess$n, res_preprocess$interpolated_x,col = "blue");
+dev.off()
+par(op)
+
+result1<- est.stats.summary(
+par = test2_res_D2_nlminb$par,
+sse = test2_res_D2_nlminb$objective,
+x = res_preprocess$interpolated_x,
+y = res_preprocess$interpolated_y)
+
+result2 <- est.summary(
+par = test2_res_D2_nlminb$par,
+sse = test2_res_D2_nlminb$objective,
+n = res_preprocess$n,
+x = res_preprocess$interpolated_x,
+y = res_preprocess$interpolated_y)
+
+
+write.table(result1,outputAfilename,row.name = F,col.name = F)
+write.table(result2,outputBfilename,row.name = F,col.name = F)
Property changes on: SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/Shellpeak.R
___________________________________________________________________
Name: svn:executable
+ *
Added: SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/Shellpeak_PK2runbyrun.R
===================================================================
--- SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/Shellpeak_PK2runbyrun.R (rev 0)
+++ SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/Shellpeak_PK2runbyrun.R 2009-04-09 21:37:19 UTC (rev 2853)
@@ -0,0 +1,64 @@
+library(stats)
+source("peak_scripts/preprocess.Enewsmooth.R")
+source("peak_scripts/peakfit_2.1.R")
+
+# get quick and dirty estimate for where baseline should be for the assessment.
+# because the cycle is 17on,10off, 63% of points should be above whatever baseline really is
+# meaning that if we sort the y-=values of the 108TR series, it should be rank #33
+
+#------------------------- the above is Uri's notes for his time series, mine below ----------------
+#--- doing PK2 run by run; there are 8 total runs, each with 226 TRs
+#--- roughly %50 are null events, rest have an event in each run
+
+allinputs <- Sys.getenv("R_SWIFT_ARGS")
+inputfilename <- noquote(strsplit(allinputs," ")[[1]][1])
+outobjectname <- noquote(strsplit(allinputs," ")[[1]][2])
+print(outobjectname)
+outputPDFfilename <- paste(outobjectname,"pdf",sep="")
+outputAfilename <- paste(outobjectname,"a",sep="")
+outputBfilename <- paste(outobjectname,"b",sep="")
+
+dd <- read.table(inputfilename)
+## make sure to change line below to fit the number of timepoints in ts
+quickbase <- (sort(dd[,1])[113])
+adjusted.input = dd[,1]-quickbase
+res_preprocess <- preprocess_1(input = adjusted.input);
+
+test2_res_D2_nlminb <- nlminb(
+ start = res_preprocess$guess,
+ objective = fcn_L2,
+ gradient = test_D_fcn_L2,
+#hessian = test_D2_fcn_L2,
+ control = list(eval.max = 50000, iter.max = 5000),
+# the lower line with fewer iterations is just for testing purposes. the upper line is the real one
+# control = list(eval.max = 100, iter.max = 50),
+ lower = rep(c(0),length(res_preprocess$guess)),
+ fcall = f,
+ n = res_preprocess$n,
+ x = res_preprocess$interpolated_x,
+ y_and_noise = res_preprocess$interpolated_y)
+
+pdf(file=outputPDFfilename)
+op <- par(mfrow=c(1,1))
+plot(res_preprocess$interpolated_x,res_preprocess$interpolated_y,xlab = "",ylab = "BOLD (1E-3)",type = "b");
+plot.est(par = test2_res_D2_nlminb$par, fcall = f, n = res_preprocess$n, res_preprocess$interpolated_x,col = "blue");
+plot.est.individual(par = test2_res_D2_nlminb$par, fcall = f, n = res_preprocess$n, res_preprocess$interpolated_x,col = "blue");
+dev.off()
+par(op)
+
+result1<- est.stats.summary(
+ par = test2_res_D2_nlminb$par,
+ sse = test2_res_D2_nlminb$objective,
+ x = res_preprocess$interpolated_x,
+ y = res_preprocess$interpolated_y)
+
+result2 <- est.summary(
+ par = test2_res_D2_nlminb$par,
+ sse = test2_res_D2_nlminb$objective,
+ n = res_preprocess$n,
+ x = res_preprocess$interpolated_x,
+ y = res_preprocess$interpolated_y)
+
+
+write.table(result1,outputAfilename,row.name = F,col.name = F)
+write.table(result2,outputBfilename,row.name = F,col.name = F)
Added: SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/peakfit_2.1.R
===================================================================
--- SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/peakfit_2.1.R (rev 0)
+++ SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/peakfit_2.1.R 2009-04-09 21:37:19 UTC (rev 2853)
@@ -0,0 +1,385 @@
+# ----------------------------------------------------------------------
+# Savitzky-Golay Algorithm
+# ----------------------------------------------------------------------
+# T2 <- sav.gol(T, fl, forder, dorder);
+#
+# Polynomial filtering method of Savitzky and Golay
+# See Numerical Recipes, 1992, Chapter 14.8, for details.
+#
+# T = vector of signals to be filtered
+# (the derivative is calculated for each ROW)
+# fl = filter length (for instance fl = 51..151)
+# forder = filter order (2 = quadratic filter, 4= quartic)
+# dorder = derivative order (0 = smoothing, 1 = first derivative, etc.)
+#
+sav.gol <- function(T, fl, forder, dorder)
+{
+ m <- length(T)
+ dorder <- dorder + 1
+
+ # -- calculate filter coefficients --
+ fc <- (fl-1)/2 # index: window left and right
+ X <- outer(-fc:fc, 0:forder, FUN="^") # polynomial terms and
+coefficients
+ Y <- pinv(X); # pseudoinverse
+
+ # -- filter via convolution and take care of the end points --
+ T2 <- convolve(T, rev(Y[dorder,]), type="o") # convolve(...)
+ T2 <- T2[(fc+1):(length(T2)-fc)]
+}
+#-----------------------------------------------------------------------
+# *** PseudoInvers of a Matrix ***
+# using singular value decomposition
+#
+pinv <- function (A)
+{
+ s <- svd(A)
+ s$v %*% diag(1/s$d) %*% t(s$u)
+}
+#-----------------------------------------------------------------------
+#
+#
+#------------gamma function used in peakfit-------------------------------------------
+#
+# | x-a1 | | (x-a1)/a2 + a3 -1 |^(a3-1)
+# y = a0*exp|- -------|*| ------------------ |
+# | a2 | | a3-1 |
+# a0 = amplitude
+# a1 = center
+# a2 = width (>0)
+# a3 = shape (> 1.01, < 167.92)
+#
+
+f_peakfit <- function(a0,a1,a2,a3,n,x) {
+ expr <- rep(c(0),length(x));
+ for (i in 1:n)
+ {
+ tmp <- rep(c(0),length(x))
+ tmp <- a0[i]*exp(-(x - a1[i])/a2[i])*((x - a1[i])/(a2[i]*(a3[i] - 1)) + 1 )^(a3[i] - 1);
+# when { b[i]/shape*(x -center[i]) + 1 } < 0, the power is meaningless. tmp[j] is "NaN" and will be replaced by zero
+ for (j in 1:length(x)) {if (tmp[j] == "NaN") {tmp[j] = 0} }
+ expr <- expr + tmp
+ }
+ eval(expr)
+}
+
+
+# define the function (a kind of simplified version from above) to be used to fit the data
+#
+# y = a*exp(-b*(x - center))*(b/shape*(x - center) + 1)^(shape)
+#
+# a = amplitude
+# b = decay parameter (correspond to 1/a2)
+# center = a1;
+# shape = a3 - 1; later shape will be fixed as 8.6, which is a common number used in most of fMRI study
+#
+#
+# -------------------------------------------------------------------------------------------
+# the following is a summation of gmma functions to fit raw data
+# assumption: linear system
+# x is the time points
+# n is the number of peaks detected by second derivatives
+# a, b and center are the parameters defined above and need to be estimated
+# a, b and center are vectors with length corresponding to the number of peaks
+#
+
+
+# shape is fixed at 8.6 in all the following analysis
+shape = 8.6
+f <- function(a,b,center,n,x) {
+ shape = 8.6;
+ expr <- rep(c(0),length(x));
+ for (i in 1:n)
+ {
+ tmp <- rep(c(0),length(x))
+ tmp <- a[i]*exp(-b[i]*(x - center[i]))*(b[i]/shape*(x -center[i]) + 1 )^shape
+# when { b[i]/shape*(x -center[i]) + 1 } < 0, the power is meaningless. tmp[j] is "NaN" and will be replaced by zero
+ for (j in 1:length(x)) {if (tmp[j] == "NaN") {tmp[j] = 0} }
+ expr <- expr + tmp
+ }
+ eval(expr)
+}
+
+# linear square fit, minmize funtion 0.5*|y - y_hat|^2.
+# Use as an input function for constraint estimation
+fcn_L2 <- function(par,n,x,y_and_noise,fcall) {
+ a <- par[1:n]
+ b <- par[(n+1):(2*n)]
+ center <- par[(2*n+1):(3*n)]
+ new_par <- c(list(a = a),list(b = b),list(center = center))
+ res <- (y_and_noise - do.call("fcall",c( as.list(new_par), list(n = n), list(x = x) ) ));
+ 0.5*t(res) %*% res
+ }
+
+
+# define the function (a kind of simplified version from above) to be used to fit the data
+#
+# y = a*exp(-b*(x - center))*(b/shape*(x - center) + 1)^(shape)
+#
+# a = amplitude
+# b = decay parameter (correspond to 1/a2)
+# center = a1;
+# shape = a3 - 1; later shape will be fixed as 8.6, which is a common number used in most of fMRI study
+#
+# df_a = exp(-b*(x - center))*(b/shape*(x - center) + 1)^(shape)
+# df_b = a*(-(x-center))*exp(-b*(x - center))*(b/shape*(x - center) + 1)^(shape) + a*exp(-b*(x - center))*shape*(b/shape*(x - center) + 1)^(shape-1)*(x - center)/shape
+# df_center = a*b*exp(-b*(x - center))*(b/shape*(x - center) + 1)^(shape) + a*exp(-b*(x - center))*shape*(b/shape*(x - center) + 1)^(shape-1)*(-b/shape)
+
+D_fcn_L2 <- function(par,n,x,y_and_noise,fcall) {
+ df_a <- matrix(rep(c(0),n*length(x)),ncol = n, nrow = length(x))
+ df_b <- matrix(rep(c(0),n*length(x)),ncol = n, nrow = length(x))
+ df_center <- matrix(rep(c(0),n*length(x)),ncol = n, nrow = length(x))
+ a <- par[1:n]
+ b <- par[(n+1):(2*n)]
+ center <- par[(2*n+1):(3*n)]
+ new_par <- c(list(a = a),list(b = b),list(center = center))
+ for (i in 1:n) {
+ df_a[,i] <- exp(-b[i]*(x - center[i]))*(b[i]/shape *(x -center[i]) + 1 )^shape
+ df_b[,i] <- a[i]*(-(x - center[i]))*exp(-b[i]*(x - center[i]))*(b[i]/shape *(x -center[i]) + 1 )^shape + a[i]*exp(-b[i]*(x - center[i]))*shape*(b[i]/shape*(x -center[i]) + 1 )^(shape - 1)*(x - center[i])/shape
+ df_center[,i] <- a[i]*b[i]*exp(-b[i]*(x - center[i]))*(b[i]/shape *(x -center[i]) + 1 )^shape + a[i]*exp(-b[i]*(x - center[i]))*shape*(b[i]/shape*(x -center[i]) + 1 )^(shape -1 )*(-b[i]/shape)
+ }
+ tmp <- cbind(df_a,df_b,df_center)
+ for (i in 1:(dim(tmp)[1]) ) {
+ for (j in 1:(dim(tmp)[2]) ) {
+ if (tmp[i,j] == "NaN") {tmp[i,j] = 0}
+ }
+ }
+ j_f <- tmp
+ y_eval <- do.call("fcall",c(list(a = a), list(b = b), list(center = center), list(n = n), list(x = x)))
+ t(j_f) %*% (y_eval - y_and_noise)
+ }
+
+#---- deriv 3-------------------------------------------------------------
+D_individual_f <- deriv3(
+ ~a*exp(-b*(x - center))*(b/shape*(x -center) + 1 )^shape,
+ c("a","b","center"),
+ function(a,b,center,x,shape) {})
+
+test_D_fcn_L2 <- function(par,n,x,y_and_noise,fcall) {
+ shape = 8.6
+ df_a <- matrix(rep(c(0),n*length(x)),ncol = n, nrow = length(x))
+ df_b <- matrix(rep(c(0),n*length(x)),ncol = n, nrow = length(x))
+ df_center <- matrix(rep(c(0),n*length(x)),ncol = n, nrow = length(x))
+ a <- par[1:n]
+ b <- par[(n+1):(2*n)]
+ center <- par[(2*n+1):(3*n)]
+ new_par <- c(list(a = a),list(b = b),list(center = center))
+ for (i in 1:n) {
+ tmp <- D_individual_f(a[i],b[i],center[i],x,shape)
+ g <- attr(tmp,"gradient")
+ df_a[,i] <- g[,1];
+ df_b[,i] <- g[,2];
+ df_center[,i] <- g[,3];
+ }
+ tmp <- cbind(df_a,df_b,df_center)
+ for (i in 1:(dim(tmp)[1]) ) {
+ for (j in 1:(dim(tmp)[2]) ) {
+ if (tmp[i,j] == "NaN") {tmp[i,j] = 0}
+ }
+ }
+ j_f <- tmp
+ y_eval <- do.call("fcall",c(list(a = a), list(b = b), list(center = center), list(n = n), list(x = x)))
+ t(j_f) %*% (y_eval - y_and_noise)
+ }
+#----------------------------------------------------------------------
+
+# 2nd derivative
+test_D2_fcn_L2 <- function(par,n,x,y_and_noise,fcall) {
+ shape = 8.6
+ df_a <- matrix(rep(c(0),n*length(x)),ncol = n, nrow = length(x))
+ df_b <- matrix(rep(c(0),n*length(x)),ncol = n, nrow = length(x))
+ df_center <- matrix(rep(c(0),n*length(x)),ncol = n, nrow = length(x))
+ a <- par[1:n]
+ b <- par[(n+1):(2*n)]
+ center <- par[(2*n+1):(3*n)]
+ new_par <- c(list(a = a),list(b = b),list(center = center))
+ f <- do.call("fcall",c(list(a = a), list(b = b), list(center = center), list(n = n), list(x = x))) - y_and_noise
+
+ for (i in 1:n) {
+ tmp <- D_individual_f(a[i],b[i],center[i],x,shape)
+ g <- attr(tmp,"gradient")
+ df_a[,i] <- g[,1];
+ df_b[,i] <- g[,2];
+ df_center[,i] <- g[,3];
+ }
+ tmp <- cbind(df_a,df_b,df_center)
+ for (i in 1:(dim(tmp)[1]) ) {
+ for (j in 1:(dim(tmp)[2]) ) {
+ if (tmp[i,j] == "NaN") {tmp[i,j] = 0}
+ }
+ }
+ j_f <- tmp
+ term1 <- t(j_f) %*% j_f
+# 1st term in the equation calculated, shall be a 3n x 3n matrix
+
+ f_D2f = matrix(rep(c(0),9*n^2),ncol = 3*n,nrow = 3*n);
+ for (i in 1:n) {
+ tmp <- D_individual_f(a[i],b[i],center[i],x,shape)
+ h <- colSums(attr(tmp,"hessian")*f,1)
+ f_D2f[i,i] = h[1,1];
+ f_D2f[n+i,i] = f_D2f[i,n+i] = h[1,2];
+ f_D2f[2*n+i,i] = f_D2f[i,2*n+i] = h[1,3];
+ f_D2f[n+i,n+i] = h[2,2];
+ f_D2f[2*n+i,n+i] = f_D2f[n+i,2*n+i] = h[2,3];
+ f_D2f[2*n+i,2*n+i] = h[3,3];
+ }
+ for (i in 1:(dim(f_D2f)[1]) ) {
+ for (j in 1:(dim(f_D2f)[2]) ) {
+ if (f_D2f[i,j] == "NaN") {f_D2f[i,j] = 0}
+ }
+ }
+
+ term1 + f_D2f
+ }
+
+#----------------------------------------------------------------------------------------
+# plot fitted curve with observed one
+plot.est <- function(par,fcall,n,x,col) {
+ a <- par[1:n]
+ b <- par[(n+1):(2*n)]
+ center <- par[(2*n+1):(3*n)]
+ new_par <- c(list(a = a),list(b = b),list(center = center))
+ y_est <- do.call("fcall",c( as.list(new_par), list(n = n), list(x = x) ))
+ lines(x,y_est,col = col,lwd = 2)
+ }
+
+# plot individual fitted curve
+plot.est.individual <- function(par,fcall,n,x,col) {
+ a <- par[1:n]
+ b <- par[(n+1):(2*n)]
+ center <- par[(2*n+1):(3*n)]
+ for (i in 1:n)
+ {
+ new_par <- c(list(a = a[i]),list(b = b[i]),list(center = center[i]))
+ y_est <- do.call("fcall",c( as.list(new_par), list(n = 1), list(x = x) ))
+ lines(x,y_est,col = col,lwd = 1)
+ }
+}
+
+
+# fit summary
+
+# fit statistical summary for constraint estimation
+est.stats.summary <- function(par,sse,x,y) {
+ n_x = length(x); # the number of observations
+ m = length(par); # the number of parameters
+ sse # sum of squares due to error
+ ave_y = mean(y); # the mean of observation
+ ssm = 0.5* t(y - ave_y) %*% (y - ave_y); # the sum of square about mean
+ r2 = 1 - sse/ssm; # R square, coefficient of determination
+ dof = n_x - m; # degree of freedom
+ r2_adjusted = 1 - sse*(n_x-1)/(ssm*(dof-1)); # adjusted R square
+ mse = sse/dof; # the mean square error
+ se = sqrt(mse); # the standard error of fit, the root mse
+ msr = (ssm - sse)/(m - 1); # mean square regression
+ f = msr/mse; # F-statistics
+ c(r2,dof,r2_adjusted,f,se)
+# c("r2" = r2,"DF" = dof,"Adj r2" = r2_adjusted,"F" = f,"standard error of fit" = se)
+ }
+
+# fit summary includes two fitting parameters (amplitude and centerfwhm and analytical area under the curve
+# -----fwhm------------
+# fwhm ~ 2.35*b^(1/2)*c
+# b & c are from gamma variate function x^b*exp(-t/c)
+# this approximation is found in ref "Event-related fMRI Contast When using Constant Interstimulus Interval: Theory and Experiment", pubished in Magnetic Resonance in Medicine 43:540-548 (2000)
+# this approximation result is a little bit different from that from peakfit
+# in our code,shape parameter is fixed at 8.6 & width related parameter is b
+# fwhm ~ 2.35*sqrt(8.6)/b
+# -----area under the curve-------
+# after some derivations, area under the curve is a/b*exp(8.6)*8.6^(-8.6)*gamma(8.6+1);
+# the derivation was confirmed with peakfit result
+est.summary_1 <- function(fcall,par,sse,n,x,y) {
+ a <- par[1:n] # estimated a
+ b <- par[(n+1):(2*n)] # estimated b
+ center <- par[(2*n+1):(3*n)] # estimated center
+# --calculate standard error of fit
+ n_x = length(x); # the number of observations
+ m = length(par); # the number of parameters
+ dof = n_x - m; # degree of freedom
+ sse # sum of squares due to error
+ mse = sse/dof; # the mean square error
+ se = sqrt(mse); # the standard error of fit, the root mse
+#---standard error of parameters, se * 1/sqrt(diag(hessian matrix of 0.5*|y-y_hat|^2))
+ std_par <- se * 1/sqrt(diag(test_D2_fcn_L2(par,n,x,y,f))) # std error of estimated parameters
+ std_a <- std_par[1:n];
+ std_b <- std_par[(n+1):(2*n)];
+ std_center <- std_par[(2*n+1):(3*n)];
+#---calculate fwhm, analytical area and percent area------
+ shape <- 8.6
+ fwhm <- matrix(c(0),ncol = 1,nrow = n);
+ area <- matrix(c(0),ncol = 1,nrow = n);
+ percent_area <- matrix(c(0),ncol = 1,nrow = n); for (i in 1:n) {
+ fwhm[i] <- 2.35*sqrt(shape)/b[i];
+ area[i] <- a[i]/b[i]*exp(shape)*shape^(-shape)*gamma(shape+1);
+ }
+ percent_area <- area/sum(area)*100;
+#---need to double check---------------------
+ std_fwhm <- std_b * 2.35*sqrt(8.6) * 1/sqrt(b^4);
+#---calculate partial F stats for each peak
+ y_est_col <- matrix(c(0),ncol = n, nrow = n_x);
+ for (i in 1:n)
+ {
+ new_par <- c(list(a = a[i]),list(b = b[i]),list(center = center[i]));
+ y_est_col[,i] <- do.call("fcall",c( as.list(new_par), list(n = 1), list(x = x) )); # each col corresponds to each peak fitted ts
+ }
+
+ p_partial <- matrix(c(0),ncol = 1, nrow = n);
+ for (i in 1:n)
+ {
+ res_wo_i <- rowSums(y_est_col[,-i]) - y; # residual without ith peak
+ ssr <- 0.5*res_wo_i %*% res_wo_i; # sum square of error of residual
+ dof_diff <- 3; # degree of freedom difference between full model and reduced model (i.e. without ith peak)
+ F_partial <- ((ssr - sse)/dof_diff)/(sse/dof); # partial F stats for ith peak
+ p_partial[i] <- 1 - pf(F_partial,dof_diff,dof);
+ }
+
+#---summary of individual results
+ summary <- data.frame(a,std_a,center,std_center,fwhm,std_fwhm,area,percent_area,p_partial);
+ name_of_row <- paste("Peak",1:n,sep="");
+ name_of_col <- c("Amplitude","e_Amplitude","Center","e_Center","FWHM","e_FWHM","Area","% Area","p_partial");
+ dimnames(summary) <- list(name_of_row,name_of_col);
+ summary
+ }
+
+est.summary <- function(par,sse,n,x,y) {
+ a <- par[1:n] # estimated a
+ b <- par[(n+1):(2*n)] # estimated b
+ center <- par[(2*n+1):(3*n)] # estimated center
+# --calculate standard error of fit
+ n_x = length(x); # the number of observations
+ m = length(par); # the number of parameters
+ dof = n_x - m; # degree of freedom
+ sse # sum of squares due to error
+ mse = sse/dof; # the mean square error
+ se = sqrt(mse); # the standard error of fit, the root mse
+#---standard error of parameters, se * 1/sqrt(diag(hessian matrix of 0.5*|y-y_hat|^2))
+ std_par <- se * 1/sqrt(diag(test_D2_fcn_L2(par,n,x,y,f))) # std error of estimated parameters
+ std_a <- std_par[1:n];
+ std_b <- std_par[(n+1):(2*n)];
+ std_center <- std_par[(2*n+1):(3*n)];
+#---calculate fwhm, analytical area and percent area------
+ shape <- 8.6
+ fwhm <- matrix(c(0),ncol = 1,nrow = n);
+ area <- matrix(c(0),ncol = 1,nrow = n);
+ percent_area <- matrix(c(0),ncol = 1,nrow = n); for (i in 1:n) {
+ fwhm[i] <- 2.35*sqrt(shape)/b[i];
+ area[i] <- a[i]/b[i]*exp(shape)*shape^(-shape)*gamma(shape+1);
+ }
+ percent_area <- area/sum(area, na.rm=T)*100;
+#---need to double check---------------------
+ std_fwhm <- std_b * 2.35*sqrt(8.6) * 1/sqrt(b^4);
+#---onset ---------------
+ onset <- center - shape/b;
+#---summary of individual results
+ summary <- data.frame(a,std_a,center,std_center,fwhm,std_fwhm,area,percent_area,onset);
+ name_of_row <- paste("Peak",1:n,sep="");
+ name_of_col <- c("Amplitude","e_Amplitude","Center","e_Center","FWHM","e_FWHM","Area","% Area","onset");
+ dimnames(summary) <- list(name_of_row,name_of_col);
+ summary
+ }
+
+
+
+
+
+
Property changes on: SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/peakfit_2.1.R
___________________________________________________________________
Name: svn:executable
+ *
Added: SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/preprocess.Enewsmooth.R
===================================================================
--- SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/preprocess.Enewsmooth.R (rev 0)
+++ SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/preprocess.Enewsmooth.R 2009-04-09 21:37:19 UTC (rev 2853)
@@ -0,0 +1,76 @@
+preprocess_1 <- function(input) {
+
+ # read file (read.table) and save as vector (as.matrix)
+ original_y <- input*1000;
+
+ # original # of points before interpolation
+ original_x <- c(1:length(original_y));
+
+ #---------calculate 2nd derivative to find peaks------------------------------------
+ #
+ # windowlength argument n
+ # p is order of polynimal
+ # D is the order of derivative
+
+ windowlength <- 3 # 3 data points covered
+ # tested 4 data points, detail info might be lost
+ order_polynomial <- 2; # includes constant, linear and quadratic terms
+ order_derivative <- 2; # need the information of second derivatives
+
+ derivative_original_y <- sav.gol(original_y, fl = windowlength, forder = order_polynomial, dorder = order_derivative);
+
+
+ # guess # of possible peaks by counting negative y derivatives and corresponding amplitude of y > noise level
+ noise_level <- 0; # so far noise level is set to 0
+
+ # guess peak location: use all negative second derivatives or use local min of second derivative
+ # peak_location = original_x[derivative_original_y < 0 & original_y > noise_level];
+
+ tmp <- rep(c(0),length(original_x));
+ for ( i in 2:(length(original_x)-1) ) {
+ if ( (derivative_original_y[i] < derivative_original_y[i-1]) & (derivative_original_y[i] < derivative_original_y[i+1]) )
+ { tmp[i] = i;}
+ }
+ peak_location = original_x[tmp > 0];
+
+ possible_peaks = length(peak_location)
+
+ #-------------------------------------------------------------------------------------
+
+
+ #-----------------interpolation data---------------------------------------------------
+ #
+ # always needs # of parameters < length of observation
+ # in the case that 2nd derivative detect many peaks, interpolation shall be used
+ stim_res = 2 # sympliy double the time points
+
+ interpolation <- ksmooth(original_x,original_y,"normal",bandwidth = 1,x.points = seq(1,length(original_x),by = 1/stim_res));
+
+ #derivative_interpolation <- sav.gol(interpolation$y, fl = windowlength, forder = order_polynomial, dorder = order_derivative);
+ # calculatae new peak location
+ new_peak_location = stim_res*(peak_location - 1) + 1;
+
+ # calculatae new peak amplitude
+
+ peak_amplitude = rep(c(0),possible_peaks);
+ for (i in 1:possible_peaks ) { peak_amplitude[i] = interpolation$y[new_peak_location[i]] }
+
+ #--------------------initializaation----------------------------------------
+ #
+ # parameter a corresponds to peak amplitude
+ a_guess <- peak_amplitude;
+
+ # parameter b corresponds to 1/peak width
+ b_guess <- rep(c(length(original_y)/possible_peaks), possible_peaks);
+
+ #real_center_guess <- interpolation$x[new_peak_location];
+ #center_guess <- b_guess*real_center_guess;
+ center_guess <- interpolation$x[new_peak_location];
+
+ guess <- c(a_guess,b_guess,center_guess);
+ n <- possible_peaks;
+
+ out <- list(guess = guess,n = n,interpolated_x = interpolation$x,interpolated_y = interpolation$y);
+ out
+}
+
\ No newline at end of file
Property changes on: SwiftApps/SIDGrid/swift/projects/andric/peakfit_pilots/PK2/peak_scripts/preprocess.Enewsmooth.R
___________________________________________________________________
Name: svn:executable
+ *
More information about the Swift-commit
mailing list