https://wiki.hevs.ch/R/api.php5?action=feedcontributions&user=Martijn.wetering&feedformat=atomR - User contributions [en]2024-03-29T13:05:29ZUser contributionsMediaWiki 1.18.1https://wiki.hevs.ch/R/index.php5/Functions/multi-exponential-model-exampleFunctions/multi-exponential-model-example2016-07-21T12:39:08Z<p>Martijn.wetering: /* multi-exponential fitting */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
==multi-exponential fitting==<br />
<br />
In this page we show an example of fitting a multi-exponential model <math>s(t) = \sum_n a_n e^{b_n t}</math>.<br />
<br />
This particular problem knows many variations and solution-methods. This is a case that focusses on the model with both b_n negative and one a_n positive and the other negative:<br />
<br />
<math>s(t) = a_d e^{-s_d t} - a_r e^{-s_r t}</math><br />
<br />
This results in a peak that is used to improve the result by adjusting the weights on both sides of the peak (the function is easily adjusted for more general cases).<br />
<br />
The used methods for exponential fitting are the integrative method and the non-linear least squares method. The first method is in fact a way to provide (nessecary) starting conditions as input for the second method. More info from:<br />
# Integrative method) http://iopscience.iop.org/article/10.1088/0957-0233/4/12/003<br />
# Review) http://scitation.aip.org/content/aip/journal/rsi/70/2/10.1063/1.1149581<br />
<br />
This example demonstrates the use of <br />
# defining functions with optional arguments<br />
# linear and non-linear models lm() and nls() from the 'stats' package.<br />
# using stop(), warning(), and break commands for debugging and <br />
<br />
Note the use of the 'plinear' algorithm in the non-linear least squares method. The model is non-linear only in the e^{b_n t} terms and not in the a_n parameters which can be found every step in the algorithm by an ordinary linear fit. The expression with cbind (creating a matrix) underlines this linear behaviour. <br />
cbind(exp(-t/tau_h),exp(-t/tau_l))<br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
<br />
#function parameters<br />
p_a <- -0.9<br />
p_b <- 1.6<br />
p_c <- 300<br />
p_d <- 24000<br />
<br />
#dataset (interval from 1 to 2000 seconds with 500 points)<br />
t <- c(1:500)*4<br />
s <- p_a*exp(-t/p_c)+p_b*exp(-t/p_d) #model<br />
sn <- s + rnorm(500)*var(s)^0.5*0.1^0.5 #model + 10% noise<br />
<br />
fit_result <- expfit(sn,t)<br />
<br />
fit_result[[4]]<br />
fit_result[[5]]<br />
<br />
</source><br />
<br />
outcome:<br />
<br />
- fitting parameters<br />
<br />
> fit_result[[4]]<br />
a_h a_l tau_h tau_l <br />
-0.9349148 1.6417528 324.5705856 17109.2731073 <br />
> fit_result[[5]]<br />
a_h_nls a_l_nls tau_h_nls tau_l_nls <br />
-0.9327765 1.6366577 320.3326143 17611.9201066 <br />
<br />
- plot<br />
<br />
[[File:Expfit_result.png]]<br />
<br />
<br />
===Function Code===<br />
<source lang='rsplus'><br />
<br />
expfit <- function(s, t , k_peak=0, control = nls.control(maxiter = 50, tol = 1e-05, minFactor = 1/1024/64^2, printEval = FALSE, warnOnly = FALSE))<br />
{<br />
<br />
#testing input<br />
if (!is.numeric(s)) {stop("signal must be a number vector")}<br />
if (!is.numeric(t)) {stop("time must be a number vector")}<br />
if (length(s) != length(t)) {stop("signal and time must be equal length vectors")}<br />
<br />
#statistics of domain (if k_peak is set manually then this step is skipped)<br />
if (k_peak == 0) {k_peak <- which.max(s)}<br />
n <- length(s)<br />
<br />
#integrating 2 times <br />
# si <- sintegral(t,s)$cdf$y-sintegral(t,s)$cdf$y[n]<br />
# sii <- sintegral(t,si)$cdf$y-sintegral(t,si)$cdf$y[n]<br />
<br />
#integration works better with trapezium rule <br />
si[1]<-0<br />
for (i in 2:n) {<br />
si[i] <- si[i-1]+0.5*(s[i]+s[i-1])*(t[i]-t[i-1])<br />
}<br />
sii[1]<-0<br />
for (i in 2:n) {<br />
sii[i] <- sii[i-1]+0.5*(si[i]+si[i-1])*(t[i]-t[i-1])<br />
}<br />
<br />
for (loop in 1:5) { <br />
#loop a few times while adjusting t_peak<br />
print(paste0(loop,": trying integrative method with peak at k = ",k_peak, " and t =",t[k_peak]))<br />
<br />
<br />
weights = rep(1,n) #starting with uniform weights<br />
if ((k_peak>20) && ((n-k_peak)>20)) { #adjusting weights according to number of points before and after the peak<br />
weights[1:k_peak] <- n/k_peak<br />
weights[-c(1:k_peak)] <- n/(n-k_peak)<br />
}<br />
<br />
# solving Tittelbach Helmrich method for finding tau<br />
mod1<-lm(s~1+t+si+sii,data=as.data.frame(list(t=t,s=s,si=si,sii=sii)),weights=weights)<br />
RSQ <- 1-sum(resid(mod1)^2)/sum((fitted(mod1)-mean(fitted(mod1)))^2)<br />
<br />
if (RSQ<0.85) {warning("low quality fit for bi-exponential model")}<br />
<br />
#tau are zero's of quadratic function 1/tau^2+qb/tau+qc=0<br />
qb=mod1$coefficients[3] <br />
qc=-mod1$coefficients[4]<br />
tau_h <- 1 / as.numeric((-qb+sqrt(qb^2-4*qc))/2)<br />
tau_l <- 1 / as.numeric((-qb-sqrt(qb^2-4*qc))/2)<br />
<br />
# solving bi-exponential fit with tau_h and tau_l<br />
mod2<-lm(s~0+ah+al,data=as.data.frame(list(t=t,s=s,ah=exp(-t/tau_h),al=exp(-t/tau_l))),weights=weights)<br />
<br />
# coefficients of fit<br />
a_h=mod2$coefficients[1]<br />
a_l=mod2$coefficients[2]<br />
<br />
if (tau_l<0) {<br />
warning("negative time scale detected. The model is not likely from two exponentials")<br />
break<br />
}<br />
<br />
#new peak value<br />
t_peak = log((-a_h/tau_h)/(a_l/tau_l))/(1/tau_h-1/tau_l)<br />
k_peak = which.min((t_peak-t)^2)<br />
}<br />
<br />
<br />
for (loop in 1:5) { <br />
#loop a few times while adjusting t_peak<br />
print(paste0(loop,": trying nls with peak at k = ",k_peak, " and t =",t[k_peak]))<br />
<br />
<br />
weights = rep(1,n) #starting with uniform weights<br />
if ((k_peak>20) && ((n-k_peak)>20)) { #adjusting weights according to number of points before and after the peak<br />
weights[1:k_peak] <- n/k_peak<br />
weights[-c(1:k_peak)] <- n/(n-k_peak)<br />
}<br />
<br />
# solving nls<br />
mod_nls <- nls(s~cbind(exp(-t/tau_h),exp(-t/tau_l)),data=as.data.frame(list(t=t,s=s)),start=list(tau_h=tau_h,tau_l=tau_l),control=control,algorithm="plinear",weights=weights)<br />
<br />
# coefficients of fit<br />
coeff <- coef(mod_nls)<br />
tau_h_nls = coeff[1]<br />
tau_l_nls = coeff[2]<br />
a_h_nls = coeff[3]<br />
a_l_nls = coeff[4]<br />
<br />
if ((tau_h_nls<0)||(tau_l_nls<0)) {<br />
warning("negative time scale detected. The model is not likely from two exponentials")<br />
break<br />
}<br />
<br />
#new peak value<br />
t_peak = log((-a_h_nls/tau_h_nls)/(a_l_nls/tau_l_nls))/(1/tau_h_nls-1/tau_l_nls)<br />
k_peak = which.min((t_peak-t)^2)<br />
}<br />
<br />
plot(t,s)<br />
lines(t,(a_h*exp(-t/tau_h)+a_l*exp(-t/tau_l)),col=2)<br />
lines(t,(a_h_nls*exp(-t/tau_h_nls)+a_l_nls*exp(-t/tau_l_nls)),col=3)<br />
<br />
<br />
<br />
modpars_a<-as.numeric(list(a_h<-a_h,a_l<-a_l,tau_h<-tau_h,tau_l<-tau_l))<br />
modpars_b<-as.numeric(list(a_h_nls<-a_h_nls,a_l_nls<-a_l_nls,tau_h_nls<-tau_h_nls,tau_l_nls<-tau_l_nls))<br />
names(modpars_a) <- c("a_h","a_l","tau_h","tau_l")<br />
names(modpars_b) <- c("a_h_nls","a_l_nls","tau_h_nls","tau_l_nls")<br />
ret <- list(mod1,mod2,mod_nls,modpars_a,modpars_b)<br />
<br />
RSQ <- 1-sum(resid(mod2)^2)/sum((fitted(mod2)-mean(fitted(mod2)))^2)<br />
print(paste0("integrate S/N estimate = ",RSQ/(1-RSQ)))<br />
RSQ <- 1-sum(resid(mod_nls)^2)/sum((fitted(mod_nls)-mean(fitted(mod_nls)))^2)<br />
print(paste0("nls S/N estimate = ",RSQ/(1-RSQ)))<br />
<br />
return(ret)<br />
}<br />
<br />
<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Functions/multi-exponential-model-exampleFunctions/multi-exponential-model-example2016-07-21T12:38:33Z<p>Martijn.wetering: /* multi-exponential fitting */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
==multi-exponential fitting==<br />
<br />
In this page we show an example of fitting a multi-exponential model <math>s(t) = \sum_n a_n e^{b_n t}</math>.<br />
<br />
This particular problem knows many variations and solution-methods. This is a case that focusses on the model with both b_n negative and one a_n positive and the other negative <br />
<math>s(t) = a_d e^{-s_d t} - a_r e^{-s_r t}</math><br />
This results in a peak that is used to improve the result by adjusting the weights on both sides of the peak (the function is easily adjusted for more general cases).<br />
<br />
The used methods for exponential fitting are the integrative method and the non-linear least squares method. The first method is in fact a way to provide (nessecary) starting conditions as input for the second method. More info from:<br />
# Integrative method) http://iopscience.iop.org/article/10.1088/0957-0233/4/12/003<br />
# Review) http://scitation.aip.org/content/aip/journal/rsi/70/2/10.1063/1.1149581<br />
<br />
This example demonstrates the use of <br />
# defining functions with optional arguments<br />
# linear and non-linear models lm() and nls() from the 'stats' package.<br />
# using stop(), warning(), and break commands for debugging and <br />
<br />
Note the use of the 'plinear' algorithm in the non-linear least squares method. The model is non-linear only in the e^{b_n t} terms and not in the a_n parameters which can be found every step in the algorithm by an ordinary linear fit. The expression with cbind (creating a matrix) underlines this linear behaviour. <br />
cbind(exp(-t/tau_h),exp(-t/tau_l))<br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
<br />
#function parameters<br />
p_a <- -0.9<br />
p_b <- 1.6<br />
p_c <- 300<br />
p_d <- 24000<br />
<br />
#dataset (interval from 1 to 2000 seconds with 500 points)<br />
t <- c(1:500)*4<br />
s <- p_a*exp(-t/p_c)+p_b*exp(-t/p_d) #model<br />
sn <- s + rnorm(500)*var(s)^0.5*0.1^0.5 #model + 10% noise<br />
<br />
fit_result <- expfit(sn,t)<br />
<br />
fit_result[[4]]<br />
fit_result[[5]]<br />
<br />
</source><br />
<br />
outcome:<br />
<br />
- fitting parameters<br />
<br />
> fit_result[[4]]<br />
a_h a_l tau_h tau_l <br />
-0.9349148 1.6417528 324.5705856 17109.2731073 <br />
> fit_result[[5]]<br />
a_h_nls a_l_nls tau_h_nls tau_l_nls <br />
-0.9327765 1.6366577 320.3326143 17611.9201066 <br />
<br />
- plot<br />
<br />
[[File:Expfit_result.png]]<br />
<br />
<br />
===Function Code===<br />
<source lang='rsplus'><br />
<br />
expfit <- function(s, t , k_peak=0, control = nls.control(maxiter = 50, tol = 1e-05, minFactor = 1/1024/64^2, printEval = FALSE, warnOnly = FALSE))<br />
{<br />
<br />
#testing input<br />
if (!is.numeric(s)) {stop("signal must be a number vector")}<br />
if (!is.numeric(t)) {stop("time must be a number vector")}<br />
if (length(s) != length(t)) {stop("signal and time must be equal length vectors")}<br />
<br />
#statistics of domain (if k_peak is set manually then this step is skipped)<br />
if (k_peak == 0) {k_peak <- which.max(s)}<br />
n <- length(s)<br />
<br />
#integrating 2 times <br />
# si <- sintegral(t,s)$cdf$y-sintegral(t,s)$cdf$y[n]<br />
# sii <- sintegral(t,si)$cdf$y-sintegral(t,si)$cdf$y[n]<br />
<br />
#integration works better with trapezium rule <br />
si[1]<-0<br />
for (i in 2:n) {<br />
si[i] <- si[i-1]+0.5*(s[i]+s[i-1])*(t[i]-t[i-1])<br />
}<br />
sii[1]<-0<br />
for (i in 2:n) {<br />
sii[i] <- sii[i-1]+0.5*(si[i]+si[i-1])*(t[i]-t[i-1])<br />
}<br />
<br />
for (loop in 1:5) { <br />
#loop a few times while adjusting t_peak<br />
print(paste0(loop,": trying integrative method with peak at k = ",k_peak, " and t =",t[k_peak]))<br />
<br />
<br />
weights = rep(1,n) #starting with uniform weights<br />
if ((k_peak>20) && ((n-k_peak)>20)) { #adjusting weights according to number of points before and after the peak<br />
weights[1:k_peak] <- n/k_peak<br />
weights[-c(1:k_peak)] <- n/(n-k_peak)<br />
}<br />
<br />
# solving Tittelbach Helmrich method for finding tau<br />
mod1<-lm(s~1+t+si+sii,data=as.data.frame(list(t=t,s=s,si=si,sii=sii)),weights=weights)<br />
RSQ <- 1-sum(resid(mod1)^2)/sum((fitted(mod1)-mean(fitted(mod1)))^2)<br />
<br />
if (RSQ<0.85) {warning("low quality fit for bi-exponential model")}<br />
<br />
#tau are zero's of quadratic function 1/tau^2+qb/tau+qc=0<br />
qb=mod1$coefficients[3] <br />
qc=-mod1$coefficients[4]<br />
tau_h <- 1 / as.numeric((-qb+sqrt(qb^2-4*qc))/2)<br />
tau_l <- 1 / as.numeric((-qb-sqrt(qb^2-4*qc))/2)<br />
<br />
# solving bi-exponential fit with tau_h and tau_l<br />
mod2<-lm(s~0+ah+al,data=as.data.frame(list(t=t,s=s,ah=exp(-t/tau_h),al=exp(-t/tau_l))),weights=weights)<br />
<br />
# coefficients of fit<br />
a_h=mod2$coefficients[1]<br />
a_l=mod2$coefficients[2]<br />
<br />
if (tau_l<0) {<br />
warning("negative time scale detected. The model is not likely from two exponentials")<br />
break<br />
}<br />
<br />
#new peak value<br />
t_peak = log((-a_h/tau_h)/(a_l/tau_l))/(1/tau_h-1/tau_l)<br />
k_peak = which.min((t_peak-t)^2)<br />
}<br />
<br />
<br />
for (loop in 1:5) { <br />
#loop a few times while adjusting t_peak<br />
print(paste0(loop,": trying nls with peak at k = ",k_peak, " and t =",t[k_peak]))<br />
<br />
<br />
weights = rep(1,n) #starting with uniform weights<br />
if ((k_peak>20) && ((n-k_peak)>20)) { #adjusting weights according to number of points before and after the peak<br />
weights[1:k_peak] <- n/k_peak<br />
weights[-c(1:k_peak)] <- n/(n-k_peak)<br />
}<br />
<br />
# solving nls<br />
mod_nls <- nls(s~cbind(exp(-t/tau_h),exp(-t/tau_l)),data=as.data.frame(list(t=t,s=s)),start=list(tau_h=tau_h,tau_l=tau_l),control=control,algorithm="plinear",weights=weights)<br />
<br />
# coefficients of fit<br />
coeff <- coef(mod_nls)<br />
tau_h_nls = coeff[1]<br />
tau_l_nls = coeff[2]<br />
a_h_nls = coeff[3]<br />
a_l_nls = coeff[4]<br />
<br />
if ((tau_h_nls<0)||(tau_l_nls<0)) {<br />
warning("negative time scale detected. The model is not likely from two exponentials")<br />
break<br />
}<br />
<br />
#new peak value<br />
t_peak = log((-a_h_nls/tau_h_nls)/(a_l_nls/tau_l_nls))/(1/tau_h_nls-1/tau_l_nls)<br />
k_peak = which.min((t_peak-t)^2)<br />
}<br />
<br />
plot(t,s)<br />
lines(t,(a_h*exp(-t/tau_h)+a_l*exp(-t/tau_l)),col=2)<br />
lines(t,(a_h_nls*exp(-t/tau_h_nls)+a_l_nls*exp(-t/tau_l_nls)),col=3)<br />
<br />
<br />
<br />
modpars_a<-as.numeric(list(a_h<-a_h,a_l<-a_l,tau_h<-tau_h,tau_l<-tau_l))<br />
modpars_b<-as.numeric(list(a_h_nls<-a_h_nls,a_l_nls<-a_l_nls,tau_h_nls<-tau_h_nls,tau_l_nls<-tau_l_nls))<br />
names(modpars_a) <- c("a_h","a_l","tau_h","tau_l")<br />
names(modpars_b) <- c("a_h_nls","a_l_nls","tau_h_nls","tau_l_nls")<br />
ret <- list(mod1,mod2,mod_nls,modpars_a,modpars_b)<br />
<br />
RSQ <- 1-sum(resid(mod2)^2)/sum((fitted(mod2)-mean(fitted(mod2)))^2)<br />
print(paste0("integrate S/N estimate = ",RSQ/(1-RSQ)))<br />
RSQ <- 1-sum(resid(mod_nls)^2)/sum((fitted(mod_nls)-mean(fitted(mod_nls)))^2)<br />
print(paste0("nls S/N estimate = ",RSQ/(1-RSQ)))<br />
<br />
return(ret)<br />
}<br />
<br />
<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Functions/multi-exponential-model-exampleFunctions/multi-exponential-model-example2016-07-21T12:16:47Z<p>Martijn.wetering: /* multi-exponential fitting */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
==multi-exponential fitting==<br />
<br />
In this page we show an example of fitting a multi-exponential model <math>s(t) = \sum_n a_n e^{b_n t}</math>.<br />
<br />
This particular problem knows many variations and solution-methods. This is a case that focusses on the model with both b_n negative and one a_n positive and the other negative (this results in a peak that is used to improve the result by adjusting the weights on both sides of the peak): <math>s(t) = a_d e^{s_d t} - a_r e^{s_r t}</math>. <br />
<br />
The used methods for exponential fitting are the integrative method and the non-linear least squares method. The first method is in fact a way to provide (nessecary) starting conditions as input for the second method. More info from:<br />
# Integrative method) http://iopscience.iop.org/article/10.1088/0957-0233/4/12/003<br />
# Review) http://scitation.aip.org/content/aip/journal/rsi/70/2/10.1063/1.1149581<br />
<br />
This example demonstrates the use of <br />
# defining functions with optional arguments<br />
# linear and non-linear models lm() and nls() from the 'stats' package.<br />
# using stop(), warning(), and break commands for debugging and <br />
<br />
Note the use of the 'plinear' algorithm in the non-linear least squares method. The model is non-linear only in the e^{b_n t} terms and not in the a_n parameters which can be found every step in the algorithm by an ordinary linear fit. The expression with cbind (creating a matrix) underlines this linear behaviour. <br />
cbind(exp(-t/tau_h),exp(-t/tau_l))<br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
<br />
#function parameters<br />
p_a <- -0.9<br />
p_b <- 1.6<br />
p_c <- 300<br />
p_d <- 24000<br />
<br />
#dataset (interval from 1 to 2000 seconds with 500 points)<br />
t <- c(1:500)*4<br />
s <- p_a*exp(-t/p_c)+p_b*exp(-t/p_d) #model<br />
sn <- s + rnorm(500)*var(s)^0.5*0.1^0.5 #model + 10% noise<br />
<br />
fit_result <- expfit(sn,t)<br />
<br />
fit_result[[4]]<br />
fit_result[[5]]<br />
<br />
</source><br />
<br />
outcome:<br />
<br />
- fitting parameters<br />
<br />
> fit_result[[4]]<br />
a_h a_l tau_h tau_l <br />
-0.9349148 1.6417528 324.5705856 17109.2731073 <br />
> fit_result[[5]]<br />
a_h_nls a_l_nls tau_h_nls tau_l_nls <br />
-0.9327765 1.6366577 320.3326143 17611.9201066 <br />
<br />
- plot<br />
<br />
[[File:Expfit_result.png]]<br />
<br />
<br />
===Function Code===<br />
<source lang='rsplus'><br />
<br />
expfit <- function(s, t , k_peak=0, control = nls.control(maxiter = 50, tol = 1e-05, minFactor = 1/1024/64^2, printEval = FALSE, warnOnly = FALSE))<br />
{<br />
<br />
#testing input<br />
if (!is.numeric(s)) {stop("signal must be a number vector")}<br />
if (!is.numeric(t)) {stop("time must be a number vector")}<br />
if (length(s) != length(t)) {stop("signal and time must be equal length vectors")}<br />
<br />
#statistics of domain (if k_peak is set manually then this step is skipped)<br />
if (k_peak == 0) {k_peak <- which.max(s)}<br />
n <- length(s)<br />
<br />
#integrating 2 times <br />
# si <- sintegral(t,s)$cdf$y-sintegral(t,s)$cdf$y[n]<br />
# sii <- sintegral(t,si)$cdf$y-sintegral(t,si)$cdf$y[n]<br />
<br />
#integration works better with trapezium rule <br />
si[1]<-0<br />
for (i in 2:n) {<br />
si[i] <- si[i-1]+0.5*(s[i]+s[i-1])*(t[i]-t[i-1])<br />
}<br />
sii[1]<-0<br />
for (i in 2:n) {<br />
sii[i] <- sii[i-1]+0.5*(si[i]+si[i-1])*(t[i]-t[i-1])<br />
}<br />
<br />
for (loop in 1:5) { <br />
#loop a few times while adjusting t_peak<br />
print(paste0(loop,": trying integrative method with peak at k = ",k_peak, " and t =",t[k_peak]))<br />
<br />
<br />
weights = rep(1,n) #starting with uniform weights<br />
if ((k_peak>20) && ((n-k_peak)>20)) { #adjusting weights according to number of points before and after the peak<br />
weights[1:k_peak] <- n/k_peak<br />
weights[-c(1:k_peak)] <- n/(n-k_peak)<br />
}<br />
<br />
# solving Tittelbach Helmrich method for finding tau<br />
mod1<-lm(s~1+t+si+sii,data=as.data.frame(list(t=t,s=s,si=si,sii=sii)),weights=weights)<br />
RSQ <- 1-sum(resid(mod1)^2)/sum((fitted(mod1)-mean(fitted(mod1)))^2)<br />
<br />
if (RSQ<0.85) {warning("low quality fit for bi-exponential model")}<br />
<br />
#tau are zero's of quadratic function 1/tau^2+qb/tau+qc=0<br />
qb=mod1$coefficients[3] <br />
qc=-mod1$coefficients[4]<br />
tau_h <- 1 / as.numeric((-qb+sqrt(qb^2-4*qc))/2)<br />
tau_l <- 1 / as.numeric((-qb-sqrt(qb^2-4*qc))/2)<br />
<br />
# solving bi-exponential fit with tau_h and tau_l<br />
mod2<-lm(s~0+ah+al,data=as.data.frame(list(t=t,s=s,ah=exp(-t/tau_h),al=exp(-t/tau_l))),weights=weights)<br />
<br />
# coefficients of fit<br />
a_h=mod2$coefficients[1]<br />
a_l=mod2$coefficients[2]<br />
<br />
if (tau_l<0) {<br />
warning("negative time scale detected. The model is not likely from two exponentials")<br />
break<br />
}<br />
<br />
#new peak value<br />
t_peak = log((-a_h/tau_h)/(a_l/tau_l))/(1/tau_h-1/tau_l)<br />
k_peak = which.min((t_peak-t)^2)<br />
}<br />
<br />
<br />
for (loop in 1:5) { <br />
#loop a few times while adjusting t_peak<br />
print(paste0(loop,": trying nls with peak at k = ",k_peak, " and t =",t[k_peak]))<br />
<br />
<br />
weights = rep(1,n) #starting with uniform weights<br />
if ((k_peak>20) && ((n-k_peak)>20)) { #adjusting weights according to number of points before and after the peak<br />
weights[1:k_peak] <- n/k_peak<br />
weights[-c(1:k_peak)] <- n/(n-k_peak)<br />
}<br />
<br />
# solving nls<br />
mod_nls <- nls(s~cbind(exp(-t/tau_h),exp(-t/tau_l)),data=as.data.frame(list(t=t,s=s)),start=list(tau_h=tau_h,tau_l=tau_l),control=control,algorithm="plinear",weights=weights)<br />
<br />
# coefficients of fit<br />
coeff <- coef(mod_nls)<br />
tau_h_nls = coeff[1]<br />
tau_l_nls = coeff[2]<br />
a_h_nls = coeff[3]<br />
a_l_nls = coeff[4]<br />
<br />
if ((tau_h_nls<0)||(tau_l_nls<0)) {<br />
warning("negative time scale detected. The model is not likely from two exponentials")<br />
break<br />
}<br />
<br />
#new peak value<br />
t_peak = log((-a_h_nls/tau_h_nls)/(a_l_nls/tau_l_nls))/(1/tau_h_nls-1/tau_l_nls)<br />
k_peak = which.min((t_peak-t)^2)<br />
}<br />
<br />
plot(t,s)<br />
lines(t,(a_h*exp(-t/tau_h)+a_l*exp(-t/tau_l)),col=2)<br />
lines(t,(a_h_nls*exp(-t/tau_h_nls)+a_l_nls*exp(-t/tau_l_nls)),col=3)<br />
<br />
<br />
<br />
modpars_a<-as.numeric(list(a_h<-a_h,a_l<-a_l,tau_h<-tau_h,tau_l<-tau_l))<br />
modpars_b<-as.numeric(list(a_h_nls<-a_h_nls,a_l_nls<-a_l_nls,tau_h_nls<-tau_h_nls,tau_l_nls<-tau_l_nls))<br />
names(modpars_a) <- c("a_h","a_l","tau_h","tau_l")<br />
names(modpars_b) <- c("a_h_nls","a_l_nls","tau_h_nls","tau_l_nls")<br />
ret <- list(mod1,mod2,mod_nls,modpars_a,modpars_b)<br />
<br />
RSQ <- 1-sum(resid(mod2)^2)/sum((fitted(mod2)-mean(fitted(mod2)))^2)<br />
print(paste0("integrate S/N estimate = ",RSQ/(1-RSQ)))<br />
RSQ <- 1-sum(resid(mod_nls)^2)/sum((fitted(mod_nls)-mean(fitted(mod_nls)))^2)<br />
print(paste0("nls S/N estimate = ",RSQ/(1-RSQ)))<br />
<br />
return(ret)<br />
}<br />
<br />
<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Functions/multi-exponential-model-exampleFunctions/multi-exponential-model-example2016-07-21T12:16:19Z<p>Martijn.wetering: /* multi-exponential fitting */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
==multi-exponential fitting==<br />
<br />
In this page we show an example of fitting a multi-exponential model <math>s(t) = \sum_n a_n e^{b_n t}</math>.<br />
<br />
This particular problem knows many variations and solution-methods. This is a case that focusses on the model with both b_n negative and one a_n positive and the other negative (this results in a peak that is used to improve the result by adjusting the weights on both sides of the peak): <math>s(t) = a_d e^{s_d t} - a_r e^{s_r t}</math>. <br />
<br />
The used methods for exponential fitting are the integrative method and the non-linear least squares method. The first method is in fact a way to provide (nessecary) starting conditions as input for the second method. More info from:<br />
# Integrative method) http://iopscience.iop.org/article/10.1088/0957-0233/4/12/003<br />
# Review) http://scitation.aip.org/content/aip/journal/rsi/70/2/10.1063/1.1149581<br />
<br />
This example demonstrates the use of <br />
# defining functions with optional arguments<br />
# linear and non-linear models lm() and nls() from the 'stats' package.<br />
# using stop(), warning(), and break commands for debugging and <br />
<br />
Note the use of the 'plinear' algorithm in the non-linear least squares method. The model is non-linear only in the e^{b_n t} terms and not in the a_n parameters which can be found every step in the algorithm by an ordinary linear fit. The expression with cbind underlines this linear behaviour. <br />
cbind(exp(-t/tau_h),exp(-t/tau_l))<br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
<br />
#function parameters<br />
p_a <- -0.9<br />
p_b <- 1.6<br />
p_c <- 300<br />
p_d <- 24000<br />
<br />
#dataset (interval from 1 to 2000 seconds with 500 points)<br />
t <- c(1:500)*4<br />
s <- p_a*exp(-t/p_c)+p_b*exp(-t/p_d) #model<br />
sn <- s + rnorm(500)*var(s)^0.5*0.1^0.5 #model + 10% noise<br />
<br />
fit_result <- expfit(sn,t)<br />
<br />
fit_result[[4]]<br />
fit_result[[5]]<br />
<br />
</source><br />
<br />
outcome:<br />
<br />
- fitting parameters<br />
<br />
> fit_result[[4]]<br />
a_h a_l tau_h tau_l <br />
-0.9349148 1.6417528 324.5705856 17109.2731073 <br />
> fit_result[[5]]<br />
a_h_nls a_l_nls tau_h_nls tau_l_nls <br />
-0.9327765 1.6366577 320.3326143 17611.9201066 <br />
<br />
- plot<br />
<br />
[[File:Expfit_result.png]]<br />
<br />
<br />
===Function Code===<br />
<source lang='rsplus'><br />
<br />
expfit <- function(s, t , k_peak=0, control = nls.control(maxiter = 50, tol = 1e-05, minFactor = 1/1024/64^2, printEval = FALSE, warnOnly = FALSE))<br />
{<br />
<br />
#testing input<br />
if (!is.numeric(s)) {stop("signal must be a number vector")}<br />
if (!is.numeric(t)) {stop("time must be a number vector")}<br />
if (length(s) != length(t)) {stop("signal and time must be equal length vectors")}<br />
<br />
#statistics of domain (if k_peak is set manually then this step is skipped)<br />
if (k_peak == 0) {k_peak <- which.max(s)}<br />
n <- length(s)<br />
<br />
#integrating 2 times <br />
# si <- sintegral(t,s)$cdf$y-sintegral(t,s)$cdf$y[n]<br />
# sii <- sintegral(t,si)$cdf$y-sintegral(t,si)$cdf$y[n]<br />
<br />
#integration works better with trapezium rule <br />
si[1]<-0<br />
for (i in 2:n) {<br />
si[i] <- si[i-1]+0.5*(s[i]+s[i-1])*(t[i]-t[i-1])<br />
}<br />
sii[1]<-0<br />
for (i in 2:n) {<br />
sii[i] <- sii[i-1]+0.5*(si[i]+si[i-1])*(t[i]-t[i-1])<br />
}<br />
<br />
for (loop in 1:5) { <br />
#loop a few times while adjusting t_peak<br />
print(paste0(loop,": trying integrative method with peak at k = ",k_peak, " and t =",t[k_peak]))<br />
<br />
<br />
weights = rep(1,n) #starting with uniform weights<br />
if ((k_peak>20) && ((n-k_peak)>20)) { #adjusting weights according to number of points before and after the peak<br />
weights[1:k_peak] <- n/k_peak<br />
weights[-c(1:k_peak)] <- n/(n-k_peak)<br />
}<br />
<br />
# solving Tittelbach Helmrich method for finding tau<br />
mod1<-lm(s~1+t+si+sii,data=as.data.frame(list(t=t,s=s,si=si,sii=sii)),weights=weights)<br />
RSQ <- 1-sum(resid(mod1)^2)/sum((fitted(mod1)-mean(fitted(mod1)))^2)<br />
<br />
if (RSQ<0.85) {warning("low quality fit for bi-exponential model")}<br />
<br />
#tau are zero's of quadratic function 1/tau^2+qb/tau+qc=0<br />
qb=mod1$coefficients[3] <br />
qc=-mod1$coefficients[4]<br />
tau_h <- 1 / as.numeric((-qb+sqrt(qb^2-4*qc))/2)<br />
tau_l <- 1 / as.numeric((-qb-sqrt(qb^2-4*qc))/2)<br />
<br />
# solving bi-exponential fit with tau_h and tau_l<br />
mod2<-lm(s~0+ah+al,data=as.data.frame(list(t=t,s=s,ah=exp(-t/tau_h),al=exp(-t/tau_l))),weights=weights)<br />
<br />
# coefficients of fit<br />
a_h=mod2$coefficients[1]<br />
a_l=mod2$coefficients[2]<br />
<br />
if (tau_l<0) {<br />
warning("negative time scale detected. The model is not likely from two exponentials")<br />
break<br />
}<br />
<br />
#new peak value<br />
t_peak = log((-a_h/tau_h)/(a_l/tau_l))/(1/tau_h-1/tau_l)<br />
k_peak = which.min((t_peak-t)^2)<br />
}<br />
<br />
<br />
for (loop in 1:5) { <br />
#loop a few times while adjusting t_peak<br />
print(paste0(loop,": trying nls with peak at k = ",k_peak, " and t =",t[k_peak]))<br />
<br />
<br />
weights = rep(1,n) #starting with uniform weights<br />
if ((k_peak>20) && ((n-k_peak)>20)) { #adjusting weights according to number of points before and after the peak<br />
weights[1:k_peak] <- n/k_peak<br />
weights[-c(1:k_peak)] <- n/(n-k_peak)<br />
}<br />
<br />
# solving nls<br />
mod_nls <- nls(s~cbind(exp(-t/tau_h),exp(-t/tau_l)),data=as.data.frame(list(t=t,s=s)),start=list(tau_h=tau_h,tau_l=tau_l),control=control,algorithm="plinear",weights=weights)<br />
<br />
# coefficients of fit<br />
coeff <- coef(mod_nls)<br />
tau_h_nls = coeff[1]<br />
tau_l_nls = coeff[2]<br />
a_h_nls = coeff[3]<br />
a_l_nls = coeff[4]<br />
<br />
if ((tau_h_nls<0)||(tau_l_nls<0)) {<br />
warning("negative time scale detected. The model is not likely from two exponentials")<br />
break<br />
}<br />
<br />
#new peak value<br />
t_peak = log((-a_h_nls/tau_h_nls)/(a_l_nls/tau_l_nls))/(1/tau_h_nls-1/tau_l_nls)<br />
k_peak = which.min((t_peak-t)^2)<br />
}<br />
<br />
plot(t,s)<br />
lines(t,(a_h*exp(-t/tau_h)+a_l*exp(-t/tau_l)),col=2)<br />
lines(t,(a_h_nls*exp(-t/tau_h_nls)+a_l_nls*exp(-t/tau_l_nls)),col=3)<br />
<br />
<br />
<br />
modpars_a<-as.numeric(list(a_h<-a_h,a_l<-a_l,tau_h<-tau_h,tau_l<-tau_l))<br />
modpars_b<-as.numeric(list(a_h_nls<-a_h_nls,a_l_nls<-a_l_nls,tau_h_nls<-tau_h_nls,tau_l_nls<-tau_l_nls))<br />
names(modpars_a) <- c("a_h","a_l","tau_h","tau_l")<br />
names(modpars_b) <- c("a_h_nls","a_l_nls","tau_h_nls","tau_l_nls")<br />
ret <- list(mod1,mod2,mod_nls,modpars_a,modpars_b)<br />
<br />
RSQ <- 1-sum(resid(mod2)^2)/sum((fitted(mod2)-mean(fitted(mod2)))^2)<br />
print(paste0("integrate S/N estimate = ",RSQ/(1-RSQ)))<br />
RSQ <- 1-sum(resid(mod_nls)^2)/sum((fitted(mod_nls)-mean(fitted(mod_nls)))^2)<br />
print(paste0("nls S/N estimate = ",RSQ/(1-RSQ)))<br />
<br />
return(ret)<br />
}<br />
<br />
<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Functions/multi-exponential-model-exampleFunctions/multi-exponential-model-example2016-07-21T12:05:06Z<p>Martijn.wetering: Created page with "{{private}} {{TOC right}} ==multi-exponential fitting== In this page we show an example of fitting a multi-exponential model <math>s(t) = \sum_n a_n e^{b_n t}</math>. This ..."</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
==multi-exponential fitting==<br />
<br />
In this page we show an example of fitting a multi-exponential model <math>s(t) = \sum_n a_n e^{b_n t}</math>.<br />
<br />
This particular problem knows many variations and solution-methods. This is a case that focusses on the model with both b_n negative and one a_n positive and the other negative (this results in a peak that is used to improve the result by adjusting the weights on both sides of the peak): <math>s(t) = a_d e^{s_d t} - a_r e^{s_r t}</math>. <br />
<br />
The used methods for exponential fitting are the integrative method and the non-linear least squares method. The first method is in fact a way to provide (nessecary) starting conditions as input for the second method. More info from:<br />
# Integrative method) http://iopscience.iop.org/article/10.1088/0957-0233/4/12/003<br />
# Review) http://scitation.aip.org/content/aip/journal/rsi/70/2/10.1063/1.1149581<br />
<br />
This example demonstrates the use of <br />
# defining functions with optional arguments<br />
# linear and non-linear models lm() and nls() from the 'stats' package.<br />
# using stop(), warning(), and break commands for debugging and <br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
<br />
#function parameters<br />
p_a <- -0.9<br />
p_b <- 1.6<br />
p_c <- 300<br />
p_d <- 24000<br />
<br />
#dataset (interval from 1 to 2000 seconds with 500 points)<br />
t <- c(1:500)*4<br />
s <- p_a*exp(-t/p_c)+p_b*exp(-t/p_d) #model<br />
sn <- s + rnorm(500)*var(s)^0.5*0.1^0.5 #model + 10% noise<br />
<br />
fit_result <- expfit(sn,t)<br />
<br />
fit_result[[4]]<br />
fit_result[[5]]<br />
<br />
</source><br />
<br />
outcome:<br />
<br />
- fitting parameters<br />
<br />
> fit_result[[4]]<br />
a_h a_l tau_h tau_l <br />
-0.9349148 1.6417528 324.5705856 17109.2731073 <br />
> fit_result[[5]]<br />
a_h_nls a_l_nls tau_h_nls tau_l_nls <br />
-0.9327765 1.6366577 320.3326143 17611.9201066 <br />
<br />
- plot<br />
<br />
[[File:Expfit_result.png]]<br />
<br />
<br />
===Function Code===<br />
<source lang='rsplus'><br />
<br />
expfit <- function(s, t , k_peak=0, control = nls.control(maxiter = 50, tol = 1e-05, minFactor = 1/1024/64^2, printEval = FALSE, warnOnly = FALSE))<br />
{<br />
<br />
#testing input<br />
if (!is.numeric(s)) {stop("signal must be a number vector")}<br />
if (!is.numeric(t)) {stop("time must be a number vector")}<br />
if (length(s) != length(t)) {stop("signal and time must be equal length vectors")}<br />
<br />
#statistics of domain (if k_peak is set manually then this step is skipped)<br />
if (k_peak == 0) {k_peak <- which.max(s)}<br />
n <- length(s)<br />
<br />
#integrating 2 times <br />
# si <- sintegral(t,s)$cdf$y-sintegral(t,s)$cdf$y[n]<br />
# sii <- sintegral(t,si)$cdf$y-sintegral(t,si)$cdf$y[n]<br />
<br />
#integration works better with trapezium rule <br />
si[1]<-0<br />
for (i in 2:n) {<br />
si[i] <- si[i-1]+0.5*(s[i]+s[i-1])*(t[i]-t[i-1])<br />
}<br />
sii[1]<-0<br />
for (i in 2:n) {<br />
sii[i] <- sii[i-1]+0.5*(si[i]+si[i-1])*(t[i]-t[i-1])<br />
}<br />
<br />
for (loop in 1:5) { <br />
#loop a few times while adjusting t_peak<br />
print(paste0(loop,": trying integrative method with peak at k = ",k_peak, " and t =",t[k_peak]))<br />
<br />
<br />
weights = rep(1,n) #starting with uniform weights<br />
if ((k_peak>20) && ((n-k_peak)>20)) { #adjusting weights according to number of points before and after the peak<br />
weights[1:k_peak] <- n/k_peak<br />
weights[-c(1:k_peak)] <- n/(n-k_peak)<br />
}<br />
<br />
# solving Tittelbach Helmrich method for finding tau<br />
mod1<-lm(s~1+t+si+sii,data=as.data.frame(list(t=t,s=s,si=si,sii=sii)),weights=weights)<br />
RSQ <- 1-sum(resid(mod1)^2)/sum((fitted(mod1)-mean(fitted(mod1)))^2)<br />
<br />
if (RSQ<0.85) {warning("low quality fit for bi-exponential model")}<br />
<br />
#tau are zero's of quadratic function 1/tau^2+qb/tau+qc=0<br />
qb=mod1$coefficients[3] <br />
qc=-mod1$coefficients[4]<br />
tau_h <- 1 / as.numeric((-qb+sqrt(qb^2-4*qc))/2)<br />
tau_l <- 1 / as.numeric((-qb-sqrt(qb^2-4*qc))/2)<br />
<br />
# solving bi-exponential fit with tau_h and tau_l<br />
mod2<-lm(s~0+ah+al,data=as.data.frame(list(t=t,s=s,ah=exp(-t/tau_h),al=exp(-t/tau_l))),weights=weights)<br />
<br />
# coefficients of fit<br />
a_h=mod2$coefficients[1]<br />
a_l=mod2$coefficients[2]<br />
<br />
if (tau_l<0) {<br />
warning("negative time scale detected. The model is not likely from two exponentials")<br />
break<br />
}<br />
<br />
#new peak value<br />
t_peak = log((-a_h/tau_h)/(a_l/tau_l))/(1/tau_h-1/tau_l)<br />
k_peak = which.min((t_peak-t)^2)<br />
}<br />
<br />
<br />
for (loop in 1:5) { <br />
#loop a few times while adjusting t_peak<br />
print(paste0(loop,": trying nls with peak at k = ",k_peak, " and t =",t[k_peak]))<br />
<br />
<br />
weights = rep(1,n) #starting with uniform weights<br />
if ((k_peak>20) && ((n-k_peak)>20)) { #adjusting weights according to number of points before and after the peak<br />
weights[1:k_peak] <- n/k_peak<br />
weights[-c(1:k_peak)] <- n/(n-k_peak)<br />
}<br />
<br />
# solving nls<br />
mod_nls <- nls(s~cbind(exp(-t/tau_h),exp(-t/tau_l)),data=as.data.frame(list(t=t,s=s)),start=list(tau_h=tau_h,tau_l=tau_l),control=control,algorithm="plinear",weights=weights)<br />
<br />
# coefficients of fit<br />
coeff <- coef(mod_nls)<br />
tau_h_nls = coeff[1]<br />
tau_l_nls = coeff[2]<br />
a_h_nls = coeff[3]<br />
a_l_nls = coeff[4]<br />
<br />
if ((tau_h_nls<0)||(tau_l_nls<0)) {<br />
warning("negative time scale detected. The model is not likely from two exponentials")<br />
break<br />
}<br />
<br />
#new peak value<br />
t_peak = log((-a_h_nls/tau_h_nls)/(a_l_nls/tau_l_nls))/(1/tau_h_nls-1/tau_l_nls)<br />
k_peak = which.min((t_peak-t)^2)<br />
}<br />
<br />
plot(t,s)<br />
lines(t,(a_h*exp(-t/tau_h)+a_l*exp(-t/tau_l)),col=2)<br />
lines(t,(a_h_nls*exp(-t/tau_h_nls)+a_l_nls*exp(-t/tau_l_nls)),col=3)<br />
<br />
<br />
<br />
modpars_a<-as.numeric(list(a_h<-a_h,a_l<-a_l,tau_h<-tau_h,tau_l<-tau_l))<br />
modpars_b<-as.numeric(list(a_h_nls<-a_h_nls,a_l_nls<-a_l_nls,tau_h_nls<-tau_h_nls,tau_l_nls<-tau_l_nls))<br />
names(modpars_a) <- c("a_h","a_l","tau_h","tau_l")<br />
names(modpars_b) <- c("a_h_nls","a_l_nls","tau_h_nls","tau_l_nls")<br />
ret <- list(mod1,mod2,mod_nls,modpars_a,modpars_b)<br />
<br />
RSQ <- 1-sum(resid(mod2)^2)/sum((fitted(mod2)-mean(fitted(mod2)))^2)<br />
print(paste0("integrate S/N estimate = ",RSQ/(1-RSQ)))<br />
RSQ <- 1-sum(resid(mod_nls)^2)/sum((fitted(mod_nls)-mean(fitted(mod_nls)))^2)<br />
print(paste0("nls S/N estimate = ",RSQ/(1-RSQ)))<br />
<br />
return(ret)<br />
}<br />
<br />
<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/File:Expfit_result.pngFile:Expfit result.png2016-07-21T11:58:36Z<p>Martijn.wetering: example of a fit to a multi-exponential model
S/N ~ 10</p>
<hr />
<div>example of a fit to a multi-exponential model<br />
S/N ~ 10</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/FunctionsFunctions2016-07-21T11:17:45Z<p>Martijn.wetering: </p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
In this section various functions are described, ranging from simplest calculations to miscellaneous functions, which cannot be classified ad hoc.<br />
<br />
= Function Sections =<br />
# [[Functions/Calculation|Calculation]]<br />
# [[Functions/Miscellaneous|Miscellaneous]]<br />
<br />
= Example Functions =<br />
# applying a [[Functions/Savitzky-Golay-filter-example|Savitzky-Golay filter]] on a 2-n matrix with x and y values as column vectors.<br />
# fitting a [[Functions/multi-exponential-model-example|multiple exponential model]] on a time series s(t).<br />
<br />
{{navNamed|left=Links|left_name=Links|up=Main_Page|up_name=Main|right=Attributes|right_name=Attributes}}<br />
<br />
[[Category:Functions]]</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Plotting/example-of-ggplot-usePlotting/example-of-ggplot-use2016-02-25T16:31:18Z<p>Martijn.wetering: /* example of use */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
<br />
An example of ggplot use. The example entails multiple concepts besides ggplot. It is a class that extends the functions of the ggplot class. It allows to personalize plots without having to re-type so much everytime. It also allows to keep the data manipulations for plotting out of sight such that the code concentrates on other manipulations (calculations, etc.)<br />
<br />
This example demonstrates the use of <br />
# ggplot<br />
# S4 classes<br />
# getting nice log scales in ggplot with the 'scales' package<br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
q <- TinusPlot() #generate<br />
<br />
<br />
# add data <br />
# note: the functions take q and output a new q (in my opinion a disadvantages in R that does not have functions that manipulate objects).<br />
<br />
q <- addLineData(q, c(1:10), c(1:10)^2, "test1")<br />
q <- addData(q, c(1:10), jitter(rep(1, 10), 10)*c(1:10)^2, "test1")<br />
<br />
q <- addLineData(q, c(1:10), c(1:10)^2.3, "test2")<br />
q <- addData(q, c(1:10), jitter(rep(1, 10), 10)*c(1:10)^2.3, "test2")<br />
<br />
# create a ggplot object (could be manipulated furter<br />
p <- plotTinus(q, log=c(TRUE, TRUE), legend=FALSE) <br />
# generate a plot<br />
p<br />
</source><br />
<br />
[[File:Gg_plot_example.png|557px|copy picture of output from above code]]<br />
<br />
===Class Code===<br />
<source lang='rsplus'><br />
require(ggplot2) # plotting<br />
require(methods) # S4 class<br />
require(scales) # logaritmic scales<br />
<br />
setClass("gg") # nessecary to use the gg S3 class in a S4 class<br />
<br />
<br />
#TinusPlot is an S4 class that parses a ggplot object based on a DATA set.<br />
# DATA a container for data that we can manipulate (add extra series) before plotting<br />
<br />
setClass("TinusPlot",<br />
slots = c(data = "data.frame"<br />
), <br />
<br />
prototype = list(<br />
)<br />
)<br />
<br />
# Upon creation of the TinusPlot class an object is created that contains <br />
# parameters that are used later with the plotTinus function to do the final plotting<br />
<br />
# internally the class defines a matrix DATA<br />
# with values for x and y <br />
# and with legend, legend2 and legendt columns<br />
# those last three columns are used to define different sets of series...<br />
<br />
# ...this is done because ggplot does not plot multiple 2xn matrixes <br />
# ggplot uses a single 2xn data matrix and additional 'factors' to define how to deal with different points<br />
<br />
<br />
TinusPlot <- function()<br />
{<br />
obj <- new("TinusPlot")<br />
<br />
#data (empty container)<br />
xp <- numeric()<br />
xl <- numeric() <br />
y_point <- numeric()<br />
y_line <- numeric()<br />
legend <- character()<br />
legend2 <- character()<br />
legendt <- character()<br />
obj@data <- as.data.frame(rbind(xp, xl, y_point, y_line, legend, legend2, legendt), stringsAsFactors = FALSE)<br />
<br />
return(obj)<br />
}<br />
<br />
# the plotTinus function is used to generate a gg class that can be plotted<br />
# several options can be defined to change the appearance<br />
# you can change this function in any way that you like to get your personal preferences<br />
<br />
setGeneric("plotTinus",<br />
function(obj, labels=c("x", "y"), size=5, textsize=12, log=c(FALSE, FALSE), theme="classic", black_dot=TRUE, legend=TRUE, diffplot=TRUE, diffline=TRUE)<br />
{<br />
standardGeneric("plotTinus")<br />
} <br />
)<br />
setMethod("plotTinus",<br />
signature("TinusPlot"),<br />
function(obj, labels=c("x", "y"), size=5, textsize=12, log=c(FALSE, FALSE), theme="classic", black_dot=TRUE, legend=TRUE, diffplot=TRUE, diffline=TRUE)<br />
{<br />
# conversion to numeric (simple cleaning up erroneous input, it does NOT detect erroneous input) <br />
obj@data$xp <- as.numeric(as.character(obj@data$xp))<br />
obj@data$xl <- as.numeric(as.character(obj@data$xl))<br />
obj@data$y_point <- as.numeric(as.character(obj@data$y_point))<br />
obj@data$y_line <- as.numeric(as.character(obj@data$y_line))<br />
<br />
#creating the ggplot object with data and labels<br />
q <- ggplot(data=obj@data) + <br />
xlab(labels[1]) +<br />
ylab(labels[2]) <br />
<br />
#adding the points layer in two different ways (same or different colour)<br />
if (diffplot) {<br />
q <- q + geom_point(aes(x=xp, y=y_point, colour = legend), size=size) <br />
} else {<br />
q <- q + geom_point(aes(x=xp, y=y_point), size=obj@size) <br />
}<br />
<br />
#adding the lines layer in two different ways (same or different colour)<br />
if (diffline) {<br />
q <- q + geom_line(aes(x=xl, y=y_line,group=legend2)) <br />
} else {<br />
q <- q + geom_line(aes(x=xl, y=y_line)) <br />
}<br />
<br />
#adding options for themes (if theme is not black or classic than the ggplot standard is used)<br />
if (black_dot) q <- q + scale_colour_grey()<br />
if (theme == "black") q <- q + theme_bw(base_size = 8)<br />
if (theme == "classic") q <- q + theme_classic(base_size=8)<br />
<br />
#textsize<br />
q <- q + theme(text = element_text(size=textsize))<br />
<br />
#adding optional choice to place legend or not<br />
if (!legend) q <- q + theme(legend.position='none')<br />
#if less than 7 curves are plotted then we plot shapes <br />
# (note ambiguity legend is a boolean inside this thread, but inside q it is a column in data)<br />
if (nlevels(obj@data$legend)<7) q <- q + aes(shape = legend)<br />
<br />
#stuff to do the logaritmic plotting<br />
if (log[1]) {<br />
ratio <- max(na.omit(log(na.omit(c(obj@data$xp,obj@data$xl)),10)))-min(na.omit(log(na.omit(c(1,obj@data$xp,obj@data$xl)),10)))<br />
order <- max(0,ceiling(ratio-1)) #reducing or increasing number of ticks to getting ticks like 10^a with a integer<br />
q <- q + scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x,n=order),labels = trans_format("log10", math_format(10^.x)))<br />
}<br />
if (log[2]) {<br />
ratio <- max(na.omit(log(na.omit(c(obj@data$y_ppoint,obj@data$y_line)),10)))-min(na.omit(log(na.omit(c(1,obj@data$y_point,obj@data$y_line)),10)))<br />
yorder <- max(0,ceiling(ratio-1)) #reducing or increasing number of ticks to getting ticks like 10^a with a integer<br />
q <- q + scale_y_log10(breaks = trans_breaks("log10", function(ya) 10^ya,n=yorder),labels = trans_format("log10", math_format(10^.x)))<br />
}<br />
#returning the result<br />
return(q + theme(legend.title=element_blank()))<br />
}<br />
)<br />
<br />
<br />
# function to expand the data set with new series<br />
<br />
# "legend" will show up in the plot legend, as well as a different colour or shape of the points<br />
# you can repeat "addData" multiple times with the same "legend", this will count as a *single* series<br />
<br />
setGeneric("addData",<br />
function(obj, xp, y_point, legend)<br />
{<br />
standardGeneric("addData")<br />
} <br />
)<br />
setMethod("addData",<br />
signature("TinusPlot"),<br />
function(obj, xp, y_point, legend)<br />
{<br />
#preparation<br />
y_line <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
xl <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
legend2 <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
legendt <- legend #legendt can be used as a column that contains info about both lines and points (without NA's)<br />
<br />
#binding the stuff<br />
adder <- cbind(xp, xl, y_point, y_line, legend, legend2, legendt)<br />
obj@data <- rbind(obj@data, adder) <br />
return(obj)<br />
}<br />
)<br />
<br />
<br />
#copy of addData but with adding data for lines instead of data for points<br />
<br />
setGeneric("addLineData",<br />
function(obj, xl, y_line, legend2)<br />
{<br />
standardGeneric("addLineData")<br />
} <br />
)<br />
setMethod("addLineData",<br />
signature("TinusPlot"),<br />
function(obj, xl, y_line, legend2)<br />
{<br />
y_point <- rep(NA,length(y_line))<br />
xp <- rep(NA,length(y_line))<br />
legend <- rep(NA,length(y_line))<br />
legendt <- legend2<br />
adder <- cbind(xp, xl, y_point, y_line, legend, legend2, legendt)<br />
obj@data <- rbind(obj@data, adder) <br />
return(obj)<br />
}<br />
)<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Plotting/example-of-ggplot-usePlotting/example-of-ggplot-use2016-02-25T16:31:04Z<p>Martijn.wetering: </p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
<br />
An example of ggplot use. The example entails multiple concepts besides ggplot. It is a class that extends the functions of the ggplot class. It allows to personalize plots without having to re-type so much everytime. It also allows to keep the data manipulations for plotting out of sight such that the code concentrates on other manipulations (calculations, etc.)<br />
<br />
This example demonstrates the use of <br />
# ggplot<br />
# S4 classes<br />
# getting nice log scales in ggplot with the 'scales' package<br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
q <- TinusPlot() #generate<br />
<br />
<br />
# add data <br />
# note: the functions take q and output a new q (in my opinion a disadvantages in R that does not have functions that manipulate objects).<br />
<br />
q <- addLineData(q, c(1:10), c(1:10)^2, "test1")<br />
q <- addData(q, c(1:10), jitter(rep(1, 10), 10)*c(1:10)^2, "test1")<br />
<br />
q <- addLineData(q, c(1:10), c(1:10)^2.3, "test2")<br />
q <- addData(q, c(1:10), jitter(rep(1, 10), 10)*c(1:10)^2.3, "test2")<br />
<br />
# create a ggplot object (could be manipulated furter<br />
p <- plotTinus(q, log=c(TRUE, TRUE), legend=FALSE) <br />
# generate a plot<br />
p<br />
</source><br />
<br />
Output:<br />
[[File:Gg_plot_example.png|557px|copy picture of output from above code]] <br />
<br />
<br />
===Class Code===<br />
<source lang='rsplus'><br />
require(ggplot2) # plotting<br />
require(methods) # S4 class<br />
require(scales) # logaritmic scales<br />
<br />
setClass("gg") # nessecary to use the gg S3 class in a S4 class<br />
<br />
<br />
#TinusPlot is an S4 class that parses a ggplot object based on a DATA set.<br />
# DATA a container for data that we can manipulate (add extra series) before plotting<br />
<br />
setClass("TinusPlot",<br />
slots = c(data = "data.frame"<br />
), <br />
<br />
prototype = list(<br />
)<br />
)<br />
<br />
# Upon creation of the TinusPlot class an object is created that contains <br />
# parameters that are used later with the plotTinus function to do the final plotting<br />
<br />
# internally the class defines a matrix DATA<br />
# with values for x and y <br />
# and with legend, legend2 and legendt columns<br />
# those last three columns are used to define different sets of series...<br />
<br />
# ...this is done because ggplot does not plot multiple 2xn matrixes <br />
# ggplot uses a single 2xn data matrix and additional 'factors' to define how to deal with different points<br />
<br />
<br />
TinusPlot <- function()<br />
{<br />
obj <- new("TinusPlot")<br />
<br />
#data (empty container)<br />
xp <- numeric()<br />
xl <- numeric() <br />
y_point <- numeric()<br />
y_line <- numeric()<br />
legend <- character()<br />
legend2 <- character()<br />
legendt <- character()<br />
obj@data <- as.data.frame(rbind(xp, xl, y_point, y_line, legend, legend2, legendt), stringsAsFactors = FALSE)<br />
<br />
return(obj)<br />
}<br />
<br />
# the plotTinus function is used to generate a gg class that can be plotted<br />
# several options can be defined to change the appearance<br />
# you can change this function in any way that you like to get your personal preferences<br />
<br />
setGeneric("plotTinus",<br />
function(obj, labels=c("x", "y"), size=5, textsize=12, log=c(FALSE, FALSE), theme="classic", black_dot=TRUE, legend=TRUE, diffplot=TRUE, diffline=TRUE)<br />
{<br />
standardGeneric("plotTinus")<br />
} <br />
)<br />
setMethod("plotTinus",<br />
signature("TinusPlot"),<br />
function(obj, labels=c("x", "y"), size=5, textsize=12, log=c(FALSE, FALSE), theme="classic", black_dot=TRUE, legend=TRUE, diffplot=TRUE, diffline=TRUE)<br />
{<br />
# conversion to numeric (simple cleaning up erroneous input, it does NOT detect erroneous input) <br />
obj@data$xp <- as.numeric(as.character(obj@data$xp))<br />
obj@data$xl <- as.numeric(as.character(obj@data$xl))<br />
obj@data$y_point <- as.numeric(as.character(obj@data$y_point))<br />
obj@data$y_line <- as.numeric(as.character(obj@data$y_line))<br />
<br />
#creating the ggplot object with data and labels<br />
q <- ggplot(data=obj@data) + <br />
xlab(labels[1]) +<br />
ylab(labels[2]) <br />
<br />
#adding the points layer in two different ways (same or different colour)<br />
if (diffplot) {<br />
q <- q + geom_point(aes(x=xp, y=y_point, colour = legend), size=size) <br />
} else {<br />
q <- q + geom_point(aes(x=xp, y=y_point), size=obj@size) <br />
}<br />
<br />
#adding the lines layer in two different ways (same or different colour)<br />
if (diffline) {<br />
q <- q + geom_line(aes(x=xl, y=y_line,group=legend2)) <br />
} else {<br />
q <- q + geom_line(aes(x=xl, y=y_line)) <br />
}<br />
<br />
#adding options for themes (if theme is not black or classic than the ggplot standard is used)<br />
if (black_dot) q <- q + scale_colour_grey()<br />
if (theme == "black") q <- q + theme_bw(base_size = 8)<br />
if (theme == "classic") q <- q + theme_classic(base_size=8)<br />
<br />
#textsize<br />
q <- q + theme(text = element_text(size=textsize))<br />
<br />
#adding optional choice to place legend or not<br />
if (!legend) q <- q + theme(legend.position='none')<br />
#if less than 7 curves are plotted then we plot shapes <br />
# (note ambiguity legend is a boolean inside this thread, but inside q it is a column in data)<br />
if (nlevels(obj@data$legend)<7) q <- q + aes(shape = legend)<br />
<br />
#stuff to do the logaritmic plotting<br />
if (log[1]) {<br />
ratio <- max(na.omit(log(na.omit(c(obj@data$xp,obj@data$xl)),10)))-min(na.omit(log(na.omit(c(1,obj@data$xp,obj@data$xl)),10)))<br />
order <- max(0,ceiling(ratio-1)) #reducing or increasing number of ticks to getting ticks like 10^a with a integer<br />
q <- q + scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x,n=order),labels = trans_format("log10", math_format(10^.x)))<br />
}<br />
if (log[2]) {<br />
ratio <- max(na.omit(log(na.omit(c(obj@data$y_ppoint,obj@data$y_line)),10)))-min(na.omit(log(na.omit(c(1,obj@data$y_point,obj@data$y_line)),10)))<br />
yorder <- max(0,ceiling(ratio-1)) #reducing or increasing number of ticks to getting ticks like 10^a with a integer<br />
q <- q + scale_y_log10(breaks = trans_breaks("log10", function(ya) 10^ya,n=yorder),labels = trans_format("log10", math_format(10^.x)))<br />
}<br />
#returning the result<br />
return(q + theme(legend.title=element_blank()))<br />
}<br />
)<br />
<br />
<br />
# function to expand the data set with new series<br />
<br />
# "legend" will show up in the plot legend, as well as a different colour or shape of the points<br />
# you can repeat "addData" multiple times with the same "legend", this will count as a *single* series<br />
<br />
setGeneric("addData",<br />
function(obj, xp, y_point, legend)<br />
{<br />
standardGeneric("addData")<br />
} <br />
)<br />
setMethod("addData",<br />
signature("TinusPlot"),<br />
function(obj, xp, y_point, legend)<br />
{<br />
#preparation<br />
y_line <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
xl <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
legend2 <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
legendt <- legend #legendt can be used as a column that contains info about both lines and points (without NA's)<br />
<br />
#binding the stuff<br />
adder <- cbind(xp, xl, y_point, y_line, legend, legend2, legendt)<br />
obj@data <- rbind(obj@data, adder) <br />
return(obj)<br />
}<br />
)<br />
<br />
<br />
#copy of addData but with adding data for lines instead of data for points<br />
<br />
setGeneric("addLineData",<br />
function(obj, xl, y_line, legend2)<br />
{<br />
standardGeneric("addLineData")<br />
} <br />
)<br />
setMethod("addLineData",<br />
signature("TinusPlot"),<br />
function(obj, xl, y_line, legend2)<br />
{<br />
y_point <- rep(NA,length(y_line))<br />
xp <- rep(NA,length(y_line))<br />
legend <- rep(NA,length(y_line))<br />
legendt <- legend2<br />
adder <- cbind(xp, xl, y_point, y_line, legend, legend2, legendt)<br />
obj@data <- rbind(obj@data, adder) <br />
return(obj)<br />
}<br />
)<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/File:Gg_plot_example.pngFile:Gg plot example.png2016-02-25T16:30:08Z<p>Martijn.wetering: </p>
<hr />
<div>example output of plotting in ggplot used in [[Plotting/example-of-ggplot-use|ggplot use]]</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/File:Gg_plot_example.pngFile:Gg plot example.png2016-02-25T16:29:34Z<p>Martijn.wetering: example output of plotting in ggplot used in
http://wiki.hevs.ch/R/index.php5/Plotting/example-of-ggplot-use
ggplot use</p>
<hr />
<div>example output of plotting in ggplot used in <br />
http://wiki.hevs.ch/R/index.php5/Plotting/example-of-ggplot-use<br />
[[Plotting/example-of-ggplot-use|ggplot use]]</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Plotting/example-of-ggplot-usePlotting/example-of-ggplot-use2016-02-25T16:26:12Z<p>Martijn.wetering: /* example of use */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
<br />
An example of ggplot use. The example entails multiple concepts besides ggplot. It is a class that extends the functions of the ggplot class. It allows to personalize plots without having to re-type so much everytime. It also allows to keep the data manipulations for plotting out of sight such that the code concentrates on other manipulations (calculations, etc.)<br />
<br />
This example demonstrates the use of <br />
# ggplot<br />
# S4 classes<br />
# getting nice log scales in ggplot with the 'scales' package<br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
q <- TinusPlot() #generate<br />
<br />
<br />
# add data <br />
# note: the functions take q and output a new q (in my opinion a disadvantages in R that does not have functions that manipulate objects).<br />
<br />
q <- addLineData(q, c(1:10), c(1:10)^2, "test1")<br />
q <- addData(q, c(1:10), jitter(rep(1, 10), 10)*c(1:10)^2, "test1")<br />
<br />
q <- addLineData(q, c(1:10), c(1:10)^2.3, "test2")<br />
q <- addData(q, c(1:10), jitter(rep(1, 10), 10)*c(1:10)^2.3, "test2")<br />
<br />
# create a ggplot object (could be manipulated furter<br />
p <- plotTinus(q, log=c(TRUE, TRUE), legend=FALSE) <br />
# generate a plot<br />
p<br />
</source><br />
<br />
===Class Code===<br />
<source lang='rsplus'><br />
require(ggplot2) # plotting<br />
require(methods) # S4 class<br />
require(scales) # logaritmic scales<br />
<br />
setClass("gg") # nessecary to use the gg S3 class in a S4 class<br />
<br />
<br />
#TinusPlot is an S4 class that parses a ggplot object based on a DATA set.<br />
# DATA a container for data that we can manipulate (add extra series) before plotting<br />
<br />
setClass("TinusPlot",<br />
slots = c(data = "data.frame"<br />
), <br />
<br />
prototype = list(<br />
)<br />
)<br />
<br />
# Upon creation of the TinusPlot class an object is created that contains <br />
# parameters that are used later with the plotTinus function to do the final plotting<br />
<br />
# internally the class defines a matrix DATA<br />
# with values for x and y <br />
# and with legend, legend2 and legendt columns<br />
# those last three columns are used to define different sets of series...<br />
<br />
# ...this is done because ggplot does not plot multiple 2xn matrixes <br />
# ggplot uses a single 2xn data matrix and additional 'factors' to define how to deal with different points<br />
<br />
<br />
TinusPlot <- function()<br />
{<br />
obj <- new("TinusPlot")<br />
<br />
#data (empty container)<br />
xp <- numeric()<br />
xl <- numeric() <br />
y_point <- numeric()<br />
y_line <- numeric()<br />
legend <- character()<br />
legend2 <- character()<br />
legendt <- character()<br />
obj@data <- as.data.frame(rbind(xp, xl, y_point, y_line, legend, legend2, legendt), stringsAsFactors = FALSE)<br />
<br />
return(obj)<br />
}<br />
<br />
# the plotTinus function is used to generate a gg class that can be plotted<br />
# several options can be defined to change the appearance<br />
# you can change this function in any way that you like to get your personal preferences<br />
<br />
setGeneric("plotTinus",<br />
function(obj, labels=c("x", "y"), size=5, textsize=12, log=c(FALSE, FALSE), theme="classic", black_dot=TRUE, legend=TRUE, diffplot=TRUE, diffline=TRUE)<br />
{<br />
standardGeneric("plotTinus")<br />
} <br />
)<br />
setMethod("plotTinus",<br />
signature("TinusPlot"),<br />
function(obj, labels=c("x", "y"), size=5, textsize=12, log=c(FALSE, FALSE), theme="classic", black_dot=TRUE, legend=TRUE, diffplot=TRUE, diffline=TRUE)<br />
{<br />
# conversion to numeric (simple cleaning up erroneous input, it does NOT detect erroneous input) <br />
obj@data$xp <- as.numeric(as.character(obj@data$xp))<br />
obj@data$xl <- as.numeric(as.character(obj@data$xl))<br />
obj@data$y_point <- as.numeric(as.character(obj@data$y_point))<br />
obj@data$y_line <- as.numeric(as.character(obj@data$y_line))<br />
<br />
#creating the ggplot object with data and labels<br />
q <- ggplot(data=obj@data) + <br />
xlab(labels[1]) +<br />
ylab(labels[2]) <br />
<br />
#adding the points layer in two different ways (same or different colour)<br />
if (diffplot) {<br />
q <- q + geom_point(aes(x=xp, y=y_point, colour = legend), size=size) <br />
} else {<br />
q <- q + geom_point(aes(x=xp, y=y_point), size=obj@size) <br />
}<br />
<br />
#adding the lines layer in two different ways (same or different colour)<br />
if (diffline) {<br />
q <- q + geom_line(aes(x=xl, y=y_line,group=legend2)) <br />
} else {<br />
q <- q + geom_line(aes(x=xl, y=y_line)) <br />
}<br />
<br />
#adding options for themes (if theme is not black or classic than the ggplot standard is used)<br />
if (black_dot) q <- q + scale_colour_grey()<br />
if (theme == "black") q <- q + theme_bw(base_size = 8)<br />
if (theme == "classic") q <- q + theme_classic(base_size=8)<br />
<br />
#textsize<br />
q <- q + theme(text = element_text(size=textsize))<br />
<br />
#adding optional choice to place legend or not<br />
if (!legend) q <- q + theme(legend.position='none')<br />
#if less than 7 curves are plotted then we plot shapes <br />
# (note ambiguity legend is a boolean inside this thread, but inside q it is a column in data)<br />
if (nlevels(obj@data$legend)<7) q <- q + aes(shape = legend)<br />
<br />
#stuff to do the logaritmic plotting<br />
if (log[1]) {<br />
ratio <- max(na.omit(log(na.omit(c(obj@data$xp,obj@data$xl)),10)))-min(na.omit(log(na.omit(c(1,obj@data$xp,obj@data$xl)),10)))<br />
order <- max(0,ceiling(ratio-1)) #reducing or increasing number of ticks to getting ticks like 10^a with a integer<br />
q <- q + scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x,n=order),labels = trans_format("log10", math_format(10^.x)))<br />
}<br />
if (log[2]) {<br />
ratio <- max(na.omit(log(na.omit(c(obj@data$y_ppoint,obj@data$y_line)),10)))-min(na.omit(log(na.omit(c(1,obj@data$y_point,obj@data$y_line)),10)))<br />
yorder <- max(0,ceiling(ratio-1)) #reducing or increasing number of ticks to getting ticks like 10^a with a integer<br />
q <- q + scale_y_log10(breaks = trans_breaks("log10", function(ya) 10^ya,n=yorder),labels = trans_format("log10", math_format(10^.x)))<br />
}<br />
#returning the result<br />
return(q + theme(legend.title=element_blank()))<br />
}<br />
)<br />
<br />
<br />
# function to expand the data set with new series<br />
<br />
# "legend" will show up in the plot legend, as well as a different colour or shape of the points<br />
# you can repeat "addData" multiple times with the same "legend", this will count as a *single* series<br />
<br />
setGeneric("addData",<br />
function(obj, xp, y_point, legend)<br />
{<br />
standardGeneric("addData")<br />
} <br />
)<br />
setMethod("addData",<br />
signature("TinusPlot"),<br />
function(obj, xp, y_point, legend)<br />
{<br />
#preparation<br />
y_line <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
xl <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
legend2 <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
legendt <- legend #legendt can be used as a column that contains info about both lines and points (without NA's)<br />
<br />
#binding the stuff<br />
adder <- cbind(xp, xl, y_point, y_line, legend, legend2, legendt)<br />
obj@data <- rbind(obj@data, adder) <br />
return(obj)<br />
}<br />
)<br />
<br />
<br />
#copy of addData but with adding data for lines instead of data for points<br />
<br />
setGeneric("addLineData",<br />
function(obj, xl, y_line, legend2)<br />
{<br />
standardGeneric("addLineData")<br />
} <br />
)<br />
setMethod("addLineData",<br />
signature("TinusPlot"),<br />
function(obj, xl, y_line, legend2)<br />
{<br />
y_point <- rep(NA,length(y_line))<br />
xp <- rep(NA,length(y_line))<br />
legend <- rep(NA,length(y_line))<br />
legendt <- legend2<br />
adder <- cbind(xp, xl, y_point, y_line, legend, legend2, legendt)<br />
obj@data <- rbind(obj@data, adder) <br />
return(obj)<br />
}<br />
)<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Plotting/example-of-ggplot-usePlotting/example-of-ggplot-use2016-02-25T16:25:54Z<p>Martijn.wetering: Created page with "{{private}} {{TOC right}} An example of ggplot use. The example entails multiple concepts besides ggplot. It is a class that extends the functions of the ggplot class. It al..."</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
<br />
An example of ggplot use. The example entails multiple concepts besides ggplot. It is a class that extends the functions of the ggplot class. It allows to personalize plots without having to re-type so much everytime. It also allows to keep the data manipulations for plotting out of sight such that the code concentrates on other manipulations (calculations, etc.)<br />
<br />
This example demonstrates the use of <br />
# ggplot<br />
# S4 classes<br />
# getting nice log scales in ggplot with the 'scales' package<br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
q <- TinusPlot() #generate<br />
<br />
<br />
# add data <br />
# note: the functions take q and output a new q (in my opinion a disadvantages in R that does not have functions that manipulate objects).<br />
<br />
q <- addLineData(q, c(1:10), c(1:10)^2, "test1")<br />
q <- addData(q, c(1:10), jitter(rep(1, 10), 10)*c(1:10)^2, "test1")<br />
<br />
q <- addLineData(q, c(1:10), c(1:10)^2.3, "test2")<br />
q <- addData(q, c(1:100), jitter(rep(1, 100), 10)*c(1:100)^2.3, "test2")<br />
<br />
# create a ggplot object (could be manipulated furter<br />
p <- plotTinus(q, log=c(TRUE, TRUE), legend=FALSE) <br />
# generate a plot<br />
p<br />
</source><br />
<br />
<br />
<br />
===Class Code===<br />
<source lang='rsplus'><br />
require(ggplot2) # plotting<br />
require(methods) # S4 class<br />
require(scales) # logaritmic scales<br />
<br />
setClass("gg") # nessecary to use the gg S3 class in a S4 class<br />
<br />
<br />
#TinusPlot is an S4 class that parses a ggplot object based on a DATA set.<br />
# DATA a container for data that we can manipulate (add extra series) before plotting<br />
<br />
setClass("TinusPlot",<br />
slots = c(data = "data.frame"<br />
), <br />
<br />
prototype = list(<br />
)<br />
)<br />
<br />
# Upon creation of the TinusPlot class an object is created that contains <br />
# parameters that are used later with the plotTinus function to do the final plotting<br />
<br />
# internally the class defines a matrix DATA<br />
# with values for x and y <br />
# and with legend, legend2 and legendt columns<br />
# those last three columns are used to define different sets of series...<br />
<br />
# ...this is done because ggplot does not plot multiple 2xn matrixes <br />
# ggplot uses a single 2xn data matrix and additional 'factors' to define how to deal with different points<br />
<br />
<br />
TinusPlot <- function()<br />
{<br />
obj <- new("TinusPlot")<br />
<br />
#data (empty container)<br />
xp <- numeric()<br />
xl <- numeric() <br />
y_point <- numeric()<br />
y_line <- numeric()<br />
legend <- character()<br />
legend2 <- character()<br />
legendt <- character()<br />
obj@data <- as.data.frame(rbind(xp, xl, y_point, y_line, legend, legend2, legendt), stringsAsFactors = FALSE)<br />
<br />
return(obj)<br />
}<br />
<br />
# the plotTinus function is used to generate a gg class that can be plotted<br />
# several options can be defined to change the appearance<br />
# you can change this function in any way that you like to get your personal preferences<br />
<br />
setGeneric("plotTinus",<br />
function(obj, labels=c("x", "y"), size=5, textsize=12, log=c(FALSE, FALSE), theme="classic", black_dot=TRUE, legend=TRUE, diffplot=TRUE, diffline=TRUE)<br />
{<br />
standardGeneric("plotTinus")<br />
} <br />
)<br />
setMethod("plotTinus",<br />
signature("TinusPlot"),<br />
function(obj, labels=c("x", "y"), size=5, textsize=12, log=c(FALSE, FALSE), theme="classic", black_dot=TRUE, legend=TRUE, diffplot=TRUE, diffline=TRUE)<br />
{<br />
# conversion to numeric (simple cleaning up erroneous input, it does NOT detect erroneous input) <br />
obj@data$xp <- as.numeric(as.character(obj@data$xp))<br />
obj@data$xl <- as.numeric(as.character(obj@data$xl))<br />
obj@data$y_point <- as.numeric(as.character(obj@data$y_point))<br />
obj@data$y_line <- as.numeric(as.character(obj@data$y_line))<br />
<br />
#creating the ggplot object with data and labels<br />
q <- ggplot(data=obj@data) + <br />
xlab(labels[1]) +<br />
ylab(labels[2]) <br />
<br />
#adding the points layer in two different ways (same or different colour)<br />
if (diffplot) {<br />
q <- q + geom_point(aes(x=xp, y=y_point, colour = legend), size=size) <br />
} else {<br />
q <- q + geom_point(aes(x=xp, y=y_point), size=obj@size) <br />
}<br />
<br />
#adding the lines layer in two different ways (same or different colour)<br />
if (diffline) {<br />
q <- q + geom_line(aes(x=xl, y=y_line,group=legend2)) <br />
} else {<br />
q <- q + geom_line(aes(x=xl, y=y_line)) <br />
}<br />
<br />
#adding options for themes (if theme is not black or classic than the ggplot standard is used)<br />
if (black_dot) q <- q + scale_colour_grey()<br />
if (theme == "black") q <- q + theme_bw(base_size = 8)<br />
if (theme == "classic") q <- q + theme_classic(base_size=8)<br />
<br />
#textsize<br />
q <- q + theme(text = element_text(size=textsize))<br />
<br />
#adding optional choice to place legend or not<br />
if (!legend) q <- q + theme(legend.position='none')<br />
#if less than 7 curves are plotted then we plot shapes <br />
# (note ambiguity legend is a boolean inside this thread, but inside q it is a column in data)<br />
if (nlevels(obj@data$legend)<7) q <- q + aes(shape = legend)<br />
<br />
#stuff to do the logaritmic plotting<br />
if (log[1]) {<br />
ratio <- max(na.omit(log(na.omit(c(obj@data$xp,obj@data$xl)),10)))-min(na.omit(log(na.omit(c(1,obj@data$xp,obj@data$xl)),10)))<br />
order <- max(0,ceiling(ratio-1)) #reducing or increasing number of ticks to getting ticks like 10^a with a integer<br />
q <- q + scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x,n=order),labels = trans_format("log10", math_format(10^.x)))<br />
}<br />
if (log[2]) {<br />
ratio <- max(na.omit(log(na.omit(c(obj@data$y_ppoint,obj@data$y_line)),10)))-min(na.omit(log(na.omit(c(1,obj@data$y_point,obj@data$y_line)),10)))<br />
yorder <- max(0,ceiling(ratio-1)) #reducing or increasing number of ticks to getting ticks like 10^a with a integer<br />
q <- q + scale_y_log10(breaks = trans_breaks("log10", function(ya) 10^ya,n=yorder),labels = trans_format("log10", math_format(10^.x)))<br />
}<br />
#returning the result<br />
return(q + theme(legend.title=element_blank()))<br />
}<br />
)<br />
<br />
<br />
# function to expand the data set with new series<br />
<br />
# "legend" will show up in the plot legend, as well as a different colour or shape of the points<br />
# you can repeat "addData" multiple times with the same "legend", this will count as a *single* series<br />
<br />
setGeneric("addData",<br />
function(obj, xp, y_point, legend)<br />
{<br />
standardGeneric("addData")<br />
} <br />
)<br />
setMethod("addData",<br />
signature("TinusPlot"),<br />
function(obj, xp, y_point, legend)<br />
{<br />
#preparation<br />
y_line <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
xl <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
legend2 <- rep(NA, length(y_point)) #no line data -> fill with NA<br />
legendt <- legend #legendt can be used as a column that contains info about both lines and points (without NA's)<br />
<br />
#binding the stuff<br />
adder <- cbind(xp, xl, y_point, y_line, legend, legend2, legendt)<br />
obj@data <- rbind(obj@data, adder) <br />
return(obj)<br />
}<br />
)<br />
<br />
<br />
#copy of addData but with adding data for lines instead of data for points<br />
<br />
setGeneric("addLineData",<br />
function(obj, xl, y_line, legend2)<br />
{<br />
standardGeneric("addLineData")<br />
} <br />
)<br />
setMethod("addLineData",<br />
signature("TinusPlot"),<br />
function(obj, xl, y_line, legend2)<br />
{<br />
y_point <- rep(NA,length(y_line))<br />
xp <- rep(NA,length(y_line))<br />
legend <- rep(NA,length(y_line))<br />
legendt <- legend2<br />
adder <- cbind(xp, xl, y_point, y_line, legend, legend2, legendt)<br />
obj@data <- rbind(obj@data, adder) <br />
return(obj)<br />
}<br />
)<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/PlottingPlotting2016-02-25T16:17:53Z<p>Martijn.wetering: </p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
The Plotting section is all about the visualization of data and the modification of graph parameters and appearances.<br />
<br />
<br />
...<br />
<br />
An example of [[Plotting/example-of-ggplot-use|ggplot use]]. The example entails multiple concepts besides ggplot. It is a class that extends the functions of the ggplot class. It allows to personalize plots without having to re-type so much everytime. It also allows to keep the data manipulations for plotting out of sight such that the code concentrates on other manipulations (calculations, etc.)</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Data/example-fitting-linear-modelsData/example-fitting-linear-models2016-02-25T16:06:24Z<p>Martijn.wetering: Created page with "{{private}} {{TOC right}} ==model selection== This is an example of an explanation returned to a question appearing in an R-course. It is not a full description of all R bu..."</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
==model selection==<br />
<br />
This is an example of an explanation returned to a question appearing in an R-course. <br />
It is not a full description of all R but demonstrates some principles. <br />
<br />
Question: "do I need to use a quadratic curve"<br />
<br />
The example shows some principles in model selection and demonstrates:<br />
# using anova to compare two models: anova(model1,model2) <br />
# explanation why anova function with single model is bad (dependency on order of terms in function)<br />
# comments on the use of t-test for selecting significant parameters (not a powerfull method if many parameters are included)<br />
# an example of the stepAIC function<br />
# an example of overfitting <br />
<br />
<source lang='rsplus'><br />
#the data<br />
BLG <- data.frame(inst=c(414.3333333,<br />
505.6666667,<br />
659.3333333,<br />
866,<br />
1049.333333,<br />
1047,<br />
1375,<br />
1251.666667,<br />
1477.666667<br />
)<br />
,blg=c(0.001,<br />
0.002,<br />
0.003,<br />
0.004,<br />
0.005,<br />
0.006,<br />
0.007,<br />
0.008,<br />
0.009<br />
)<br />
)<br />
<br />
####### not nessecary just alternative ##########<br />
# it is all basically linear algebra <br />
# so I like to write out explicitly the Matrix instead of using I() as done below<br />
data <- list(y<-BLG$inst, x<-BLG$blg, x2<-BLG$blg^2)<br />
model_alternative <- lm(y~1+x+x2, data=data)<br />
summary(model_alternative)<br />
####### not nessecary just alternative ##########<br />
<br />
####### four different models<br />
model_minimal <- lm(inst~1, data<-BLG )<br />
model_small1 <- lm(inst~1+blg, data<-BLG )<br />
model_small2 <- lm(inst~1+I(blg^2), data<-BLG )<br />
model_big <- lm(inst~1+blg+I(blg^2), data<-BLG )<br />
<br />
####### compare the change in RSS and see if it is significant with F test<br />
# F-test is comparing the ratio of RSS1 and RSS2<br />
# A small p(F>what you have) means that you have a small chance for such F or larger<br />
# thus the model with small RSS (causing big F) is significant improvement<br />
# A big p(F>what you have) means that you have a big chance for such F or larger<br />
# thus the model with small RSS (causing big F) is not a significant improvement (any random non meaningfull model could have done the same)<br />
<br />
anova(model_minimal,model_small1)<br />
# big F, reduction of RSS, with high significance (doesn't happen by chance)<br />
# thus y=a*1+b*blg is better than y=a*1<br />
anova(model_minimal,model_small2)<br />
# y=a*1+c*blg^2 is also better than y=a*1 but not in the same amount<br />
anova(model_small1,model_big)<br />
# y=a*1+b*blg+c*blg^2 is not better than y=a*1+b*blg<br />
<br />
# it is risky to do anova function with a single model<br />
# effectively it does multiple dual anova tests in order that you place the terms in the function<br />
model_big_diff <- lm(inst~1+I(blg^2)+blg, data<-BLG )<br />
anova(model_big)<br />
anova(model_big_diff)<br />
<br />
# you can also do a t-test based on the estimates of the parameters<br />
summary(model_big)<br />
# but beware that large models increase the error in the estimates (d freedom)<br />
# so it is very strict and easily rejects a value, it is not always powerfull<br />
# see for example the change in t-test when the model is smaller (and Std. Error gets smaller)<br />
summary(model_small1)<br />
<br />
<br />
<br />
# a better way than the single anova <br />
# (explained earlier, the stuff about order of terms)<br />
# is to use stepAIC<br />
# it adds (AND REMOVES) to find a model with low RSS that is not to low in df<br />
require(MASS)<br />
stepAIC(model_big,test="F",scope="~1+blg+I(blg^2)")<br />
<br />
# it is an algoritm you can start from anypoint <br />
# this time you end the same but note that the found solution is not always the same<br />
# and depends on the starting point<br />
<br />
# also, this is an explorative method.<br />
# You'd have to test the model again with new data. <br />
# The p-values are about the steps, not about the entire test<br />
# if you just feed enough terms the chance is high that something comes out of it<br />
# so you'd have to consider that when evaluating the p-value<br />
<br />
stepAIC(model_minimal,test="F",scope="~1+blg+I(blg^2)")<br />
<br />
<br />
# this anova tells you if going from linear to quadratic is good <br />
anova(model_small1,model_big)<br />
# it is not only improvement of df*RSS by a factor 1.0357<br />
<br />
# plotting<br />
plot(BLG$blg, BLG$inst)<br />
lines(BLG$blg, predict(model_small1))<br />
<br />
<br />
<br />
# this down here explains intuitively why<br />
# increasing models just to improve R^2 and get RSS to 0 is bad<br />
# if you go far enough you get 1<br />
model_extreme <- lm(inst~blg+I(blg^2)+I(blg^3)+I(blg^4)+I(blg^5)+I(blg^6)+I(blg^7)+I(blg^8), data<-BLG )<br />
anova(model_small1,model_extreme)<br />
<br />
# example of plotting with newdata for the predict function<br />
# also demonstration of the absurdity of that polynomial function<br />
plot(BLG$blg, BLG$inst)<br />
new <- data.frame(blg=c(1:100)*0.00009)<br />
lines(new$blg, predict(model_extreme,new))<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Data/FittingData/Fitting2016-02-25T15:57:08Z<p>Martijn.wetering: Created page with "{{private}} {{TOC right}} In this section... = examples = # Using Anova and lm function example of model selection, quadratic curve or ..."</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
In this section...<br />
<br />
= examples =<br />
# [[Data/example-fitting-linear-models|Using Anova and lm function]] example of model selection, quadratic curve or linear curve?<br />
<br />
{{navNamed|left=Links|left_name=Links|up=Main_Page|up_name=Main|right=Attributes|right_name=Attributes}}<br />
<br />
[[Category:Data]]</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Functions/Savitzky-Golay-filter-exampleFunctions/Savitzky-Golay-filter-example2015-12-30T11:39:11Z<p>Martijn.wetering: /* savitsky golay filter */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
<br />
==savitsky golay filter==<br />
<br />
This type of filter applies a piecewise fit to obtain a filtered value or a derivative (e.g compare with piecewise averaging of two or more points). <br />
<br />
This example demonstrates the use of <br />
# defining functions with optional arguments<br />
# foreach loop from the 'foreach' package<br />
# linear model lm from the 'stats' package.<br />
<br />
An already existing implementation of this filter exists in R, 'sgolayfilt' in the 'signal' package. A practical limitation of this function is that the length of the output vector is one less than the input vector, also a time series must be used with equally spaced time intervals. (note: for equally spaced time intervals the sgolayfilt, which uses a convolution instead of solving a linear model, is *much* faster)<br />
<br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
<br />
# ideal data<br />
x <- c(1:200) * 3.142/50<br />
y <- sin(x)<br />
<br />
# added noise<br />
y_j <- jitter(y, amount=0.5) <br />
<br />
# filtered<br />
y_f <- tinusGolay(x, y_j)<br />
<br />
# derivative<br />
y_fd <- tinusGolay(x, y_j, n=31, d=TRUE)<br />
<br />
<br />
# visualization<br />
plot(x, y_j, col=1)<br />
points(x, y, col=2)<br />
points(x, y_f, col=3)<br />
points(x, y_fd, col=4)<br />
<br />
</source><br />
<br />
[[File:Savitzky_golay_example.png|557px|copy picture of output from above code]] <br />
<br />
red: original sine function<br />
<br />
black: sine function after adding noise<br />
<br />
green: after filtering <br />
<br />
blue: derivative (note we used a larger set of data, 31 points, to average/filter the data, compared to the filter used to get the green curve)<br />
<br />
<br />
===Function Code===<br />
<source lang='rsplus'><br />
# tinusGolay<br />
# a function that takes two vectors of equal length and applies savitsky golay filter<br />
# this filter fits a linear model y=a+bx (it could be generalized to polynomial functions) <br />
#<br />
# this fitting is done in a piecewise manner (with pieces of length n)<br />
#<br />
# the optional parameter d is used to determine the derivative instead of the filtered value<br />
#<br />
# the optional parameter log is used to first calculate the logarithm and then back to the original value<br />
# this is usefull to find derivatives of functions whose logarithm behaves more smoothly <br />
# (such occurs in measurements that track several orders of magnitude)<br />
<br />
# note: this function does not yet control for false input! e.g. if x and y are different length<br />
<br />
<br />
require(foreach) # for generating the loop that does the piecewice fits <br />
<br />
<br />
tinusGolay <- function(x, y, n=11, d=FALSE)<br />
{<br />
l <- length(x)<br />
yn <- y #creating another array of same length, this will be the output<br />
<br />
#dataframe used in the fitting<br />
data<-as.data.frame(list(x=x, y=y))<br />
<br />
# the if statement differentiates between the filtered y and the derivative of y... <br />
# ...by the use of 'predict(model)' and 'model$coefficients[2]' <br />
if (d==FALSE) {<br />
#the first n/2 points are treated the same<br />
model <- lm(y~x, data=data[c(1:n),])<br />
yn[1:((n+1)/2)] <-predict(model)[1:((n+1)/2)]<br />
<br />
#the last n/2 points are treated the same<br />
model <- lm(y~x, data=data[c((l-n+1):l),])<br />
yn[(l-((n+1)/2)+1):l] <-predict(model)[((n+1)/2):n] <br />
<br />
#the middle points are treated differently by scanning from n/2 to l-n/2 and repeatingly applying a linear fit<br />
yn[c((((n+1)/2)+1):(l-(n+1)/2))] <- foreach (i = 2:(l-n), .combine='c') %do% {<br />
model <- lm(y~x,data=data[c(i:(i+n-1)),])<br />
predict(model)[[((n+1)/2)]]<br />
}<br />
}<br />
else {<br />
#the first n/2 points are treated the same<br />
model <- lm(y~x, data=data[c(1:n),])<br />
yn[1:((n+1)/2)] <-(model$coefficients[2])<br />
<br />
#the last n/2 points are treated the same<br />
model <- lm(y~x, data=data[c((l-n+1):l),])<br />
yn[(l-((n+1)/2)+1):l] <-(model$coefficients[2])<br />
<br />
#the middle points are treaded differentle by scanning from n/2 to l-n/2 and repeatingly applying a linear fit<br />
yn[c((((n+1)/2)+1):(l-(n+1)/2))] <- foreach (i = 2:(l-n), .combine='c') %do% {<br />
model <- lm(y~x, data=data[c(i:(i+n-1)),])<br />
(model$coefficients[2])<br />
} <br />
}<br />
<br />
return(yn) <br />
}<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Talk:Functions/Salitzky-Golay-filter-exampleTalk:Functions/Salitzky-Golay-filter-example2015-12-30T11:37:36Z<p>Martijn.wetering: Created page with "page had a spelling error in the name and can be removed"</p>
<hr />
<div>page had a spelling error in the name and can be removed</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Functions/Savitzky-Golay-filter-exampleFunctions/Savitzky-Golay-filter-example2015-12-30T11:36:11Z<p>Martijn.wetering: Created page with "{{private}} {{TOC right}} ==savitsky golay filter== This type of filter applies a piecewise fit to obtain a filtered value or a derivative (e.g compare with piecewise avera..."</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
<br />
==savitsky golay filter==<br />
<br />
This type of filter applies a piecewise fit to obtain a filtered value or a derivative (e.g compare with piecewise averaging of two or more points). <br />
<br />
This example demonstrates the use of <br />
# defining functions with optional arguments<br />
# foreach loop from the 'foreach' package<br />
# linear model lm from the 'stats' package.<br />
<br />
An already existing implementation of this filter exists in R, 'sgolayfilt' in the 'signal' package. A practical limitation of this function is that the length of the output vector is one less than the input vector, also a time series must be used with equally spaced time intervals.<br />
<br />
<br />
===example of use===<br />
<br />
<source lang='rsplus'><br />
<br />
# ideal data<br />
x <- c(1:200) * 3.142/50<br />
y <- sin(x)<br />
<br />
# added noise<br />
y_j <- jitter(y, amount=0.5) <br />
<br />
# filtered<br />
y_f <- tinusGolay(x, y_j)<br />
<br />
# derivative<br />
y_fd <- tinusGolay(x, y_j, n=31, d=TRUE)<br />
<br />
<br />
# visualization<br />
plot(x, y_j, col=1)<br />
points(x, y, col=2)<br />
points(x, y_f, col=3)<br />
points(x, y_fd, col=4)<br />
<br />
</source><br />
<br />
[[File:Savitzky_golay_example.png|557px|copy picture of output from above code]] <br />
<br />
red: original sine function<br />
<br />
black: sine function after adding noise<br />
<br />
green: after filtering <br />
<br />
blue: derivative (note we used a larger set of data, 31 points, to average/filter the data, compared to the filter used to get the green curve)<br />
<br />
<br />
===Function Code===<br />
<source lang='rsplus'><br />
# tinusGolay<br />
# a function that takes two vectors of equal length and applies savitsky golay filter<br />
# this filter fits a linear model y=a+bx (it could be generalized to polynomial functions) <br />
#<br />
# this fitting is done in a piecewise manner (with pieces of length n)<br />
#<br />
# the optional parameter d is used to determine the derivative instead of the filtered value<br />
#<br />
# the optional parameter log is used to first calculate the logarithm and then back to the original value<br />
# this is usefull to find derivatives of functions whose logarithm behaves more smoothly <br />
# (such occurs in measurements that track several orders of magnitude)<br />
<br />
# note: this function does not yet control for false input! e.g. if x and y are different length<br />
<br />
<br />
require(foreach) # for generating the loop that does the piecewice fits <br />
<br />
<br />
tinusGolay <- function(x, y, n=11, d=FALSE)<br />
{<br />
l <- length(x)<br />
yn <- y #creating another array of same length, this will be the output<br />
<br />
#dataframe used in the fitting<br />
data<-as.data.frame(list(x=x, y=y))<br />
<br />
# the if statement differentiates between the filtered y and the derivative of y... <br />
# ...by the use of 'predict(model)' and 'model$coefficients[2]' <br />
if (d==FALSE) {<br />
#the first n/2 points are treated the same<br />
model <- lm(y~x, data=data[c(1:n),])<br />
yn[1:((n+1)/2)] <-predict(model)[1:((n+1)/2)]<br />
<br />
#the last n/2 points are treated the same<br />
model <- lm(y~x, data=data[c((l-n+1):l),])<br />
yn[(l-((n+1)/2)+1):l] <-predict(model)[((n+1)/2):n] <br />
<br />
#the middle points are treated differently by scanning from n/2 to l-n/2 and repeatingly applying a linear fit<br />
yn[c((((n+1)/2)+1):(l-(n+1)/2))] <- foreach (i = 2:(l-n), .combine='c') %do% {<br />
model <- lm(y~x,data=data[c(i:(i+n-1)),])<br />
predict(model)[[((n+1)/2)]]<br />
}<br />
}<br />
else {<br />
#the first n/2 points are treated the same<br />
model <- lm(y~x, data=data[c(1:n),])<br />
yn[1:((n+1)/2)] <-(model$coefficients[2])<br />
<br />
#the last n/2 points are treated the same<br />
model <- lm(y~x, data=data[c((l-n+1):l),])<br />
yn[(l-((n+1)/2)+1):l] <-(model$coefficients[2])<br />
<br />
#the middle points are treaded differentle by scanning from n/2 to l-n/2 and repeatingly applying a linear fit<br />
yn[c((((n+1)/2)+1):(l-(n+1)/2))] <- foreach (i = 2:(l-n), .combine='c') %do% {<br />
model <- lm(y~x, data=data[c(i:(i+n-1)),])<br />
(model$coefficients[2])<br />
} <br />
}<br />
<br />
return(yn) <br />
}<br />
</source></div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/FunctionsFunctions2015-12-30T11:35:54Z<p>Martijn.wetering: /* Function Sections */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
In this section various functions are described, ranging from simplest calculations to miscellaneous functions, which cannot be classified ad hoc.<br />
<br />
= Function Sections =<br />
# [[Functions/Calculation|Calculation]]<br />
# [[Functions/Miscellaneous|Miscellaneous]]<br />
<br />
example functions<br />
# applying a [[Functions/Savitzky-Golay-filter-example|Savitzky-Golay filter]] on a 2-n matrix with x and y values as column vectors<br />
<br />
<br />
{{navNamed|left=Links|left_name=Links|up=Main_Page|up_name=Main|right=Attributes|right_name=Attributes}}<br />
<br />
[[Category:Functions]]</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/File:Savitzky_golay_example.pngFile:Savitzky golay example.png2015-12-29T18:44:17Z<p>Martijn.wetering: </p>
<hr />
<div>This is an example output used on the page <br />
http://wiki.hevs.ch/R/index.php5/Functions/Salitzky-Golay-filter-example<br />
<br />
# red: original sine function<br />
# black: sine function after adding noise<br />
# green: after filtering <br />
# blue: derivative</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/File:Savitzky_golay_example.pngFile:Savitzky golay example.png2015-12-29T18:43:49Z<p>Martijn.wetering: This is an example output used on the page
http://wiki.hevs.ch/R/index.php5/Functions/Salitzky-Golay-filter-example
red: original sine function
black: sine function after adding noise
green: after filtering
blue: derivative</p>
<hr />
<div>This is an example output used on the page <br />
http://wiki.hevs.ch/R/index.php5/Functions/Salitzky-Golay-filter-example<br />
<br />
red: original sine function<br />
black: sine function after adding noise<br />
green: after filtering <br />
blue: derivative</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/FunctionsFunctions2015-12-29T17:31:24Z<p>Martijn.wetering: /* Function Sections */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
In this section various functions are described, ranging from simplest calculations to miscellaneous functions, which cannot be classified ad hoc.<br />
<br />
= Function Sections =<br />
# [[Functions/Calculation|Calculation]]<br />
# [[Functions/Miscellaneous|Miscellaneous]]<br />
<br />
example functions<br />
# applying a [[Functions/Salitzky-Golay-filter-example|Salitzky-Golay filter]] on a 2-n matrix with x and y values as column vectors<br />
<br />
<br />
{{navNamed|left=Links|left_name=Links|up=Main_Page|up_name=Main|right=Attributes|right_name=Attributes}}<br />
<br />
[[Category:Functions]]</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/FunctionsFunctions2015-12-29T17:31:09Z<p>Martijn.wetering: /* Function Sections */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
In this section various functions are described, ranging from simplest calculations to miscellaneous functions, which cannot be classified ad hoc.<br />
<br />
= Function Sections =<br />
# [[Functions/Calculation|Calculation]]<br />
# [[Functions/Miscellaneous|Miscellaneous]]<br />
<br />
example functions<br />
# applying a [Functions/Salitzky-Golay-filter-example|Salitzky-Golay filter]] on a 2-n matrix with x and y values as column vectors<br />
<br />
<br />
{{navNamed|left=Links|left_name=Links|up=Main_Page|up_name=Main|right=Attributes|right_name=Attributes}}<br />
<br />
[[Category:Functions]]</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/FunctionsFunctions2015-12-29T17:29:44Z<p>Martijn.wetering: /* Function Sections */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
In this section various functions are described, ranging from simplest calculations to miscellaneous functions, which cannot be classified ad hoc.<br />
<br />
= Function Sections =<br />
# [[Functions/Calculation|Calculation]]<br />
# [[Functions/Miscellaneous|Miscellaneous]]<br />
<br />
an example function applying a Salitzky-Golay filter on a 2-n matrix with x and y values as column vectors<br />
# [[Functions/Salitzky-Golay-filter-example|Example of Salitzky-Golay filter]]<br />
<br />
<br />
{{navNamed|left=Links|left_name=Links|up=Main_Page|up_name=Main|right=Attributes|right_name=Attributes}}<br />
<br />
[[Category:Functions]]</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Talk:LinksTalk:Links2015-12-07T11:59:40Z<p>Martijn.wetering: Created page with "Maybe add https://google.github.io/styleguide/Rguide.xml and use the style as well on the pages in this wiki?"</p>
<hr />
<div>Maybe add https://google.github.io/styleguide/Rguide.xml and use the style as well on the pages in this wiki?</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Data/InvestigationData/Investigation2015-12-07T11:03:21Z<p>Martijn.wetering: /* Examples */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
<br />
==grep(''pattern,x,value=T'')==<br />
looks for a ''pattern'' in a vector ''x'' and returns patztern element or index of this element<br />
===Arguments===<br />
''pattern'' to look for, vector of investigation ''x'', ''value'' to determine index or element return<br />
===Examples===<br />
'''>''' x<-c(3,2,1,'hello')<br />
<br />
'''>''' grep(3,x,value=TRUE)<br />
<br />
'''[1]''' "3"<br />
<br />
'''>''' grep(3,x,value=FALSE)<br />
<br />
'''[1]''' 1<br />
<br />
'''>''' grep('hello',x,value=TRUE)<br />
<br />
'''[1]''' "hello"<br />
<br />
'''>''' grep('hello',x,value=FALSE)<br />
<br />
'''[1]''' 4<br />
<br />
<br />
'''>''' y<-c(3,'hello',2,1,'hello') #if the pattern occurs more often then a vector will be returned<br />
<br />
'''>''' grep('hello',y,value=FALSE)<br />
<br />
'''[1]''' 2 5<br />
<br />
==which.min(''x'') ; which.max(''x'')==<br />
provides the first location (as index) of the minimum / maximum in the vector<br />
===Arguments===<br />
numeric vector ''x''<br />
===Examples===<br />
'''>''' x<-c(3,2,1,3)<br />
<br />
'''>''' which.min(x)<br />
<br />
'''[1]''' 3<br />
<br />
'''>''' which.max(x)<br />
<br />
'''[1]''' 1<br />
<br />
==is.na(''x'')==<br />
checks objects or objects in a vector, whether they are ''NA'' and returns ''TRUE'' or ''FALSE'' <br />
===Arguments===<br />
to be tested object ''x''<br />
===Examples===<br />
'''>''' x<-c(1,2,NA,3,NA)<br />
<br />
'''>''' is.na(x)<br />
<br />
'''[1]''' FALSE FALSE TRUE FALSE TRUE<br />
<br />
==''x''==''y'' ; ''x''<''y'' ; ''x''>''y''==<br />
compares vector values at corresponding indices (equal, smaller as, greater as) and returns ''TRUE'' or ''FALSE'' for each index<br />
===Arguments===<br />
to be compared vectors ''x'' and ''y''<br />
===Examples===<br />
'''>''' x<-c(1,2,3,4)<br />
<br />
'''>''' y<-c(1,2,2,4)<br />
<br />
'''>''' x==y<br />
<br />
'''[1]''' TRUE TRUE FALSE TRUE</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Data/InvestigationData/Investigation2015-12-07T11:00:30Z<p>Martijn.wetering: /* Examples */</p>
<hr />
<div>{{private}}<br />
{{TOC right}}<br />
<br />
<br />
==grep(''pattern,x,value=T'')==<br />
looks for a ''pattern'' in a vector ''x'' and returns patztern element or index of this element<br />
===Arguments===<br />
''pattern'' to look for, vector of investigation ''x'', ''value'' to determine index or element return<br />
===Examples===<br />
'''>''' x<-c(3,2,1,'hello')<br />
<br />
'''>''' grep(3,x,value=TRUE)<br />
<br />
'''[1]''' "3"<br />
<br />
'''>''' grep(3,x,value=FALSE)<br />
<br />
'''[1]''' 1<br />
<br />
'''>''' grep('hello',x,value=TRUE)<br />
<br />
'''[1]''' "hello"<br />
<br />
'''>''' grep('hello',x,value=FALSE)<br />
<br />
'''[1]''' 4<br />
<br />
<br />
'''>''' y<-c(3,'hello',2,1,'hello') #if the pattern occurs more often then a vector will be returned<br />
<br />
'''>''' grep(3,y,value=FALSE)<br />
<br />
'''[1]''' 2 5<br />
<br />
==which.min(''x'') ; which.max(''x'')==<br />
provides the first location (as index) of the minimum / maximum in the vector<br />
===Arguments===<br />
numeric vector ''x''<br />
===Examples===<br />
'''>''' x<-c(3,2,1,3)<br />
<br />
'''>''' which.min(x)<br />
<br />
'''[1]''' 3<br />
<br />
'''>''' which.max(x)<br />
<br />
'''[1]''' 1<br />
<br />
==is.na(''x'')==<br />
checks objects or objects in a vector, whether they are ''NA'' and returns ''TRUE'' or ''FALSE'' <br />
===Arguments===<br />
to be tested object ''x''<br />
===Examples===<br />
'''>''' x<-c(1,2,NA,3,NA)<br />
<br />
'''>''' is.na(x)<br />
<br />
'''[1]''' FALSE FALSE TRUE FALSE TRUE<br />
<br />
==''x''==''y'' ; ''x''<''y'' ; ''x''>''y''==<br />
compares vector values at corresponding indices (equal, smaller as, greater as) and returns ''TRUE'' or ''FALSE'' for each index<br />
===Arguments===<br />
to be compared vectors ''x'' and ''y''<br />
===Examples===<br />
'''>''' x<-c(1,2,3,4)<br />
<br />
'''>''' y<-c(1,2,2,4)<br />
<br />
'''>''' x==y<br />
<br />
'''[1]''' TRUE TRUE FALSE TRUE</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/Main_PageMain Page2015-12-07T10:50:03Z<p>Martijn.wetering: /* Members */</p>
<hr />
<div>{{TOC_right}}<br />
<br />
This Wiki is a subset of the [https://www.r-project.org/ R programming language and environment for statistical computing and graphics].<br />
[[File:Hei_logo.png|left|200px|baseline|HESSO Valais Wallis Logo|link=http://www.hevs.ch]]<br />
It has its source in the R interest group, which was established in October 2015 at the University of Applied Sciences Western Switzerland Valais Wallis in Sion.<br />
The interest of the group lies in the communication of R related problems and solutions in chosen topics. Hence it is not the target to conclude another summary, user guide or manual, but rather to have a progressing documentation of tricks, hints and solution methods including the entire application scope of R. Given that this documentation originates from the R interest group communication, it should mirror the solutions of the most frequently appearing issues and obstacles during R usage.<br />
<br />
For the basic understanding of the R language several useful links are provided.<br />
<br />
[[File:Icontexto-inside-rss.png|link=Help:http://wiki.hevs.ch/R/index.php?title=Special:RecentChanges&feed=rss|48px]] [http://wiki.hevs.ch/R/index.php?title=Special:RecentChanges&feed=rss Get informed about all changes to this wiki by signing up to this RSS feed.]<br />
<br />
= Getting started =<br />
To make the experience on this wiki as pleasant as possible for everyone, please follow the<br />
* [[Help:Guidelines|Guidelines]].<br />
<br />
Then use the ''Navigation'' to the left to enter the different sections or use any of the following links:<br />
* [[Help:Contents|Help]]<br />
* [//www.mediawiki.org/wiki/Manual:FAQ MediaWiki FAQ]<br />
<br />
== Members ==<br />
* [[User:Helene.strese|Helene Strese]]<br />
* [[User:Frederic.truffer|Frederic Truffer]]<br />
* [[User:Martial.geiser|Martial Geiser]]<br />
* [[User:Oliver.gubler|Oliver Gubler]]<br />
* [[User:Djano.kandaswa|Djano Kandaswamy]]<br />
* [[User:Jessen.page|Jessen Page]]<br />
* [[User:Jean.decaix|Jean Decaix]]<br />
* [[User:Paulandr.salamin|Paul-Andre Salamin]]<br />
* [[User:Martijn.wetering|Martijn Weterings]]<br />
* [[User:Jeanmanu.segura|Jean-Manuel Segura]]<br />
* [[User:Stephani.karmann|Stephanie Karmann]]<br />
* [[User:Miriam.barque|Miriam Barque]]<br />
* [[User:Rene.schumann|Rene Schumann]]<br />
* [[User:Michael.barry|Michael Barry]]<br />
* [[User:Elianesa.maalouf|Eliane Samih Maalouf]]<br />
* [[User:Poveda.geovanny|Poveda Geovanny]]</div>Martijn.weteringhttps://wiki.hevs.ch/R/index.php5/User:Martijn.weteringUser:Martijn.wetering2015-12-07T10:47:20Z<p>Martijn.wetering: Created page with "'''Martijn Weterings''' I work in the food technology pilot hall of Batiment F. My current topic is ''the study or aroma release'' and it's relation with mass transport. I am..."</p>
<hr />
<div>'''Martijn Weterings'''<br />
<br />
I work in the food technology pilot hall of Batiment F. My current topic is ''the study or aroma release'' and it's relation with mass transport. I am using R now for all kinds data analysis and preparing plots for latex. An interessting point between R and my current work is that aroma research generates more data than you would think (there are 1000s of aroma's). Currently I am using lot's of csv files to store and manage a database. If anyone has a nice solution I am very happy to hear about it. In the past I was a student of ''food technology'' and ''astronomy'' where I used R in calculations and also, just for fun.<br />
<br />
My R interrests are mostly in the algebra and mathematical, matlab alternative, realm: ''computational mathematics, linear and non-linear modelling, factor analysis, exploratory data analysis''. <br />
My favorite plotting library is ''ggplot'' and a layer of my own written around it. I am a fan of the ''S4 classes'' approach, although I am still a rookie in this field so maybe I don't know better.<br />
<br />
I'd love to hear if any of the people from systems industriel know a way to read ''LabView'' tdm and tdx files with R-code or another programm that can be run from within R.</div>Martijn.wetering