notesofdabbler avatar notesofdabbler committed 2932ce4

directory structure change done

Comments (0)

Files changed (15)

01_learnR_pt1/20120804_CompApp.Rmd

+# Learning R using a Chemical Reaction Engineering Book
+
+[Chemical Reactor Analysis and Design Fundamentals](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/) by J.B. Rawlings and J. G. Ekerdt is a textbook for studying Chemical Reaction Engineering. The popular open source package [Octave](http://www.gnu.org/software/octave/) has its origins to the reaction engineering course offered by Prof. Rawlings. This book is accompanied by Octave and Matlab code for solving typical problems encountered in Reaction Engineering. 
+
+I figured that maybe one way to learn R is so see whether I can code some of th examples from this book in R. I am by no means suggesting that R can replace MATLAB/Octave for engineering problems but merely it is a way for me to learn the language.
+
+I started with the [computational appendix](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/web-appendix.pdf) listed in the book's website and am trying to work through some of the examples there. It will be good to refer to the computational appendix to follow the R code below.
+
+### Setting up a stoichiometric matrix, reaction rate vector and determining the rank
+
+```{r}
+# stoichiometric matrix
+stoi=matrix(c(0,1,0,-1,-1,1,
+       -1,1,1,-1,0,0,
+       1,0,-1,0,-1,1),
+       ncol=6,byrow=T)
+stoi
+
+# rank of the stoichiometrix matrix
+rank=qr(stoi)$rank
+rank
+
+# reaction rate vector
+r=c(1,2,3)
+```
+Given the reaction rate r=(r1,r2)', the rate of change of species concentration R is given by
+$$latex R=\nu^Tr $$
+where $\nu$ is the stoichiometrix matrix
+```{r}
+# rate of change of components
+R=t(stoi)%*%r
+R
+```
+### Example A1: Estimating reaction rates
+
+The stoichometrix matrix $\nu$ is input below
+```{r}
+# stoichiometry
+stoi=matrix(c(0,1,0,-1,-1,1,
+              -1,1,1,-1,0,0),nrow=2,byrow=T)
+stoi
+# number of species and number of reactions
+nspec=ncol(stoi)
+nr=nrow(stoi)
+
+nspec
+nr
+
+# true rxn rates
+r=c(1,2)
+r
+
+# true component rates
+R=t(stoi)%*%r
+R
+```
+
+### simulate 2000 measured component rates
+
+Add random noise (normally distributed with mean 0 and standard deviation 0.05) to true species rate vector R
+$$ R^m = R + \epsilon, \epsilon \sim N(0,0.0025) $$
+
+```{r}
+## simulate 2000 noise estimates
+e=matrix(0.05*rnorm(2000*nspec,0,1),nrow=2000,byrow=T)
+Rmeas=matrix(rep(R,2000),ncol=nspec,byrow=T)+e
+```
+The least squares estimate of reaction rate vector $\cap{r}$ is
+$$ \hat{r}=(\nu\nu^T)^{-1}{\nu}R $$
+
+```{r}
+## estimate reaction rates
+rest=solve(stoi%*%t(stoi),stoi%*%t(Rmeas))
+```
+I was trying different plot features in R and applying to this data of estimated rates
+
+I found the following function that plots scatterplot with marginal histograms
+```{r}
+# plotting scatterplot with histogram
+# downloaded from web
+#  http://www.r-bloggers.com/example-8-41-scatterplot-with-marginal-histograms/
+#
+scatterhist = function(x, y, xlab="", ylab=""){
+  par.default <- par(no.readonly=TRUE)
+  zones=matrix(c(2,0,1,3), ncol=2, byrow=TRUE)
+  layout(zones, widths=c(4/5,1/5), heights=c(1/5,4/5))
+  xhist = hist(x, plot=FALSE)
+  yhist = hist(y, plot=FALSE)
+  top = max(c(xhist$counts, yhist$counts))
+  par(mar=c(3,3,1,1))
+  plot(x,y)
+  par(mar=c(0,3,1,1))
+  barplot(xhist$counts, axes=FALSE, ylim=c(0, top), space=0)
+  par(mar=c(3,0,1,1))
+  barplot(yhist$counts, axes=FALSE, xlim=c(0, top), space=0, horiz=TRUE)
+  par(oma=c(3,3,0,0))
+  mtext(xlab, side=1, line=1, outer=TRUE, adj=0, 
+        at=0.9 * (mean(x) - min(x))/(max(x)-min(x)))
+  mtext(ylab, side=2, line=1, outer=TRUE, adj=0, 
+        at=(0.9 * (mean(y) - min(y))/(max(y) - min(y))))
+  par=par(par.default)
+}
+
+# scatter plot of reaction rates with marginal histograms
+scatterhist(t(rest)[,1],t(rest)[,2],xlab="r_1",ylab="r_2")
+```
+There is library cars that has a command 'scatterplot' for plotting scatterplot with box plots and has several other options
+In the plot below 50% and 90% ellipses are overlaid on the data
+
+```{r}
+# scatter plot of reaction rates with marginal box plots
+library(car)
+scatterplot(t(rest)[,1],t(rest)[,2],reg.line=FALSE,smooth=FALSE,ellipse=TRUE)
+```
+I tried ggplot2 also for the same plot
+
+```{r}
+# 2d contours/density with ggplot2 for reaction rates
+
+# create dataframe of reaction rates
+rest_df=data.frame(t(rest))
+names(rest_df)=c("r1","r2")
+
+library(ggplot2)
+ggplot(data=rest_df,aes(x=r1,y=r2))+geom_point()
+ggplot(data=rest_df,aes(x=r1,y=r2))+stat_density2d()
+ggplot(data=rest_df,aes(x=r1,y=r2))+stat_density2d(aes(fill=..level..),geom="polygon")
+```
+
+

02_learnR_pt2_3/20120811_CompApp_pt_2_3.Rmd

+# Learning R using a Chemical Reaction Engineering Book (Part 2)
+
+In case you missed part 1, you can view it [here](http://rdabbler.blogspot.com/2012/08/learning-r-using-chemical-reaction.html). In this part, I tried to recreate the examples in sections A.2.2 and A.2.3 of the [computational appendix](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/web-appendix.pdf) in the [reaction engineering book](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/) (by Rawlings and Ekerdt). 
+
+## Solving a nonlinear system of equations
+
+This example involves determining reaction equilibrium conditions by solving the following system of nonlinear equations.
+
+$$latex
+ \begin{aligned}
+   PK_1y_Iy_B-y_{P1}=0, \\
+   PK_2y_Iy_B-y_{P2}=0
+ \end{aligned}   
+$$
+
+The relation between the variables $latex y_I,y_B,y_{P1},y_{P2}$ and extent of reactions $latex x_1,x_2$ are:
+$$latex
+  \begin{aligned}
+    y_I=\frac{y_{I0}-x_1-x_2}{1-x_1-x_2}, \\
+    y_B=\frac{y_{B0}-x_1-x_2}{1-x_1-x_2}, \\ 
+    y_{P1}=\frac{y_{p10}+x_1}{1-x_1-x_2}, \\
+    y_{P2}=\frac{y_{p20}+x_2}{1-x_1-x_2}
+  \end{aligned}
+$$  
+
+Here I have used R package [rootSolve](http://cran.r-project.org/web/packages/rootSolve/index.html) for solving the above set of equations. The library is loaded and the functions to be solved are defined in the R function fns.
+
+```{r}
+# load library rootSolve
+library(rootSolve)
+
+# function defining F(x)=0
+fns=function(x){
+  K1=108; K2=284; P=2.5
+  yI0=0.5; yB0=0.5; yP10=0; yP20=0;
+  d=1-x[1]-x[2]
+  yI=(yI0-x[1]-x[2])/d
+  yB=(yB0-x[1]-x[2])/d
+  yP1=(yP10+x[1])/d
+  yP2=(yP20+x[2])/d
+  F1=P*K1*yI*yB-yP1
+  F2=P*K2*yI*yB-yP2
+  c(F1=F1,F2=F2)
+}
+```
+Next, an initial guess of (0.2,0.2) is set for the variables and the equations are solved using the function multiroot (from package rootSolve)
+
+```{r}
+# initial guess for x
+xinit=c(0.2,0.2)
+
+# solve the equations
+xans=multiroot(f=fns,start=xinit)
+
+# object returned by multiroot
+xans
+
+# solution to the equations
+xans$root
+```
+The solution to the equations is accessed from the variable xans$root which in this case is (0.1334,0.3507)
+
+## Function Minimization
+
+Here the function to be minimized is
+
+
+$$latex
+  G=-(x_1lnK_1+x_2lnK_2)+(1-x_1-x_2)lnP+(y_{I0}-x_1-x_2)ln(y_I)+(y_{B0}-x_1-x_2)ln(y_B)+(y_{P10}+x_1)ln(y_{p1})+(y_{P20}+x_2)ln(y_{P2})
+$$
+
+The following constraints need to be satisfied for $latex x_1$ and $latex x_2$
+
+$$latex
+  0 \leq x_1 \leq 0.5, 0 \leq x_2 \leq 0.5, x_1+x_2 \leq 0.5
+$$latex
+
+First I used [constrOptim](http://stat.ethz.ch/R-manual/R-patched/library/stats/html/constrOptim.html) function for minimization. We need to specify the function to be minimized
+
+```{r}
+# function to be minimized
+eval_f0=function(x){
+  dg1=-3.72e3; dg2=-4.49e3; T=400; R=1.987; P=2.5
+  K1=exp(-dg1/(R*T)); K2=exp(-dg2/(R*T))
+  
+  yI0=0.5; yB0=0.5; yP10=0; yP20=0;
+  d=1-x[1]-x[2]
+  yI=(yI0-x[1]-x[2])/d
+  yB=(yB0-x[1]-x[2])/d
+  yP1=(yP10+x[1])/d
+  yP2=(yP20+x[2])/d
+  
+  f=-(x[1]*log(K1)+x[2]*log(K2))+(1-x[1]-x[2])*log(P)+yI*d*log(yI)+
+       yB*d*log(yB)+yP1*d*log(yP1)+yP2*d*log(yP2)
+  
+  return(f)
+}
+```
+The constraints need to be specified in the form $latex Ax-b \geq 0$
+
+```{r}
+#  constraint
+A=matrix(c(-1,-1,1,0,-1,0,0,1,0,-1),ncol=2,byrow=TRUE)
+A
+b=c(-0.5,0,-0.5,0,-0.5)
+b
+```
+Next, the function is minimized using constrOptim (starting from an initial guess of (0.2,0.2)). Here Nelder-Mead method is used since BFGS method requires specifying the gradient of the function by the user. R taskview [optimization](http://cran.r-project.org/web/views/Optimization.html) lists other options.
+
+```{r}
+# initial guess
+xinit=c(0.2,0.2)
+
+# minimize function subject to bounds and constraints
+xans2=constrOptim(theta=xinit, f=eval_f0, grad=NULL, ui=A, ci=b,
+            method = "Nelder-Mead")
+xans2
+```
+The solution can be accessed from xans2$par and is (0.1331,0.3509)
+
+Next, I also tried function minimization with ConstrOptim.nl from the package [alabama](http://cran.r-project.org/web/packages/alabama/index.html). Here the constraints are specified in terms of $latex h(x) \geq 0$. Even if gradient is not supplied, this function will estimate it using finite-differences.
+
+Definition of constraints in the format for constrOptim.nl
+
+```{r}
+# load library alabama
+library(alabama)
+library(numDeriv)
+
+h_ineq=function(x){
+   h=rep(NA,1)
+   h[1]=-x[1]-x[2]+0.5
+   h[2]=x[1]
+   h[3]=-x[1]+0.5
+   h[4]=x[2]
+   h[5]=-x[2]+0.5
+   return(h)
+}
+
+xans3=constrOptim.nl(par=xinit,fn=eval_f0,hin=h_ineq)
+
+xans3
+```
+The solution can be accessed from xans3$par and is (0.1332,0.3508)
+
+Since this is a 2 dimensional problem, the solution can also be visualized using a contour plot of the function.
+
+```{r message=FALSE,warning=FALSE}
+# Region of interest: 0.01<=x1<=0.49, 0.01<=x2<=0.49
+x1=seq(0.01,0.49,by=0.01)
+x2=seq(0.01,0.49,by=0.01)
+
+# vectorizing function eval_f0 so that it can be evaluted in the outer function
+fcont=function(x,y) eval_f0(c(x,y))
+fcontv=Vectorize(fcont,SIMPLIFY=FALSE)
+z=outer(x1,x2,fcontv)
+
+# filled.coutour and contour are overlaid with the minimum point (0.133,0.351)
+filled.contour(x1,x2,z,xlab="x1",ylab="x2",
+               plot.axes={axis(1); axis(2); contour(x1,x2,z,levels=c(-3,-2.5,-2,-1.5,-1,0),vfont=c("sans serif","bold"),labcex=1,lwd=2,add=TRUE); 
+                          points(0.133,0.351,pch=15,col="blue")})
+```
+
+MATLAB/Octave functions for solving nonlinear equations (fsolve and fmincon) have been used in Chemical Engineering computations for a long time and are robust. R has traditionally not been used in this domain. So it is hard to say how the functions I have used in this blog will perform across the range of problems encountered in Reaction Engineering.
+
+This post was created with R 2.15.1 and knitr 0.7.1 (converted markdown to html)

03_learnR_pt4/20121209_CompApp_pt4.Rmd

+# Learning R using a Chemical Reaction Engineering Book (Part 3)
+
+In case you missed the previous parts, the links to them are listed below:
+* [Part1](http://rdabbler.blogspot.com/2012/08/learning-r-using-chemical-reaction.html). 
+* [Part2](http://rdabbler.blogspot.com/2012/08/learning-r-using-chemical-reaction_11.html)
+
+In this part, I tried to recreate the examples in sections A.3.1 of the [computational appendix](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/web-appendix.pdf) in the [reaction engineering book](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/) (by Rawlings and Ekerdt). 
+
+## Solving a system of ordinary differential equations
+
+This example involves reaction (Benzene pyrolysis) in a plug flow reactor. The actual reactions happening are:
+
+$$ (Rxn1) \; 2B = D + H $$
+$$ (Rxn2) \; B + D = T+H $$
+
+
+The rate of each reaction is given by:
+$$ r_1=k_1(c_B^2-\frac{c_Dc_H}{K_1}) $$
+$$ r_2=k_2(c_Bc_D-\frac{c_Tc_H}{K_2}) $$
+
+The feed to the reactor consists of 60kmol/hr of Benzene (B). The temperature of the reactor is $T=1033K$ and the pressure is $P=1atm$.The rate constants and equilibrium constants for this example are:
+$$ k_1=7\times 10^5\; L/mol.hr,\; k_2=4\times 10^5 \; L/mol.hr,\; K_1=0.31,\; K_2=0.48 $$
+
+```{r, tidy=FALSE}
+# set working directory
+setwd('~/R/RxnEngg_Rawlings')
+
+# load libraries
+library(deSolve)
+
+
+# Appendix A.3.1: Solution of Differential Equations
+
+# Benzene pyrolysis example
+
+# Parameters
+# NBf - feed benzene rate - mol/h
+# R - Universal gas constant
+# T - Reactor temperature K
+# P - Reactor pressure atm
+# k1 - rxn1 forward rate constant L/mol.h
+# k2 - rxn2 forward rate constant L/mol.h
+# Keq1 - rxn1 equilibrium constant
+# Keq2 - rxn2 equilibrium constant
+pars=c(
+NBf=60e3,  
+R=0.08205,  
+T=1033, 
+P=1, 
+k1=7e5, 
+k2=4e5, 
+Keq1=0.31, 
+Keq2=0.48  
+)
+```
+
+The governing equations for conversion versus volume in a plug flow reactor is based on extent of each of the reactions:
+$$ \frac{d\epsilon_1}{dV}=r_1,\; \frac{d\epsilon_2}{dV}=r_2 $$
+The initial conditions (corresponding to feed conditions $N_B(0)=60kmol/h,\;N_D(0)=N_H(0)=N_T(0)=0$) are that the extent of reaction is zero. 
+$$ \epsilon_1(0)=0, \; \epsilon_2(0)=0 $$
+
+The flow rates of each component along the reactor volume can be calculated from reaction extent
+$$ N_B=N_B(0)-2\epsilon_1-\epsilon_2, \; N_D=\epsilon_1-\epsilon_2, \; N_H=\epsilon_1+\epsilon_2, \; N_T=\epsilon_2 $$ 
+
+These are setup in a function that can be passed to an ODE solver. In this case the ODE solver we use lsode from the R package [deSolve](http://desolve.r-forge.r-project.org/). The inputs to the function are:
+* Variable over which the integration is done (Volume in this case)
+* The state variables of the system (Extent of the two reactions)
+* Parameters that are needed for description of the system (Rate constants, Temperature, Pressure, etc.)
+The output from this function is the rate of change as described by the equations previously.
+
+```{r, tidy=FALSE}
+# function that will be passed to odesolver
+# vol is the variable over which the system is integrated (equivalent of time in batch reactions)
+# ext is the extent of reactions 1 and 2
+# params are the parameters passed to the system
+rxnrate=function(vol,ext,params) {
+     with(as.list(c(ext,params)),{
+        NB=NBf-2*ext1-ext2
+        ND=ext1-ext2
+        NH=ext1+ext2
+        NT=ext2
+        Q=NBf*R*T/P
+        cB=NB/Q
+        cD=ND/Q
+        cT=NT/Q
+        cH=NH/Q
+        dext1=k1*(cB*cB-cD*cH/Keq1)
+        dext2=k2*(cB*cD-cT*cH/Keq2)
+        return(list(c(dext1=dext1,dext2=dext2)))
+     })
+}
+```
+
+Since the reaction start only after the feed enters the reactor, the extent of reaction is zero for both reactions at the beginning of the reactor (V=0L). The set of volumes where the concentration and reaction extent is computed is chosen in this case to be from 0L to 1600L at every 50L. The ODE solver lsode from [deSolve](http://desolve.r-forge.r-project.org/) package is used to solve the system of equations.
+
+```{r, tidy=FALSE}
+# initial extent of reaction (zero in this case for both reactions)
+extinit=c(ext1=0,ext2=0)
+# Volumes where the concentration is reported (in this case 0 to 1600L at every 50L)
+vols=seq(0,1600,length=50)
+# Solution of the set of differential equations using lsode solver in deSolve package
+extout=lsode(times=vols,y=extinit,func=rxnrate,parms=pars)
+```
+**extout** contains the extent of reaction vs volume data. That is used to compute mole fraction and conversion at different volumes along the reactor.
+
+```{r, tidy=FALSE}
+# Calcuation of mole fraction and conversion from extent of reaction at different volumes
+extoutdf=data.frame(extout)
+NBf=pars["NBf"]
+extoutdf$conv=(extoutdf$ext1*2+extoutdf$ext2)/NBf
+extoutdf$yB=(NBf-2*extoutdf$ext1-extoutdf$ext2)/NBf
+extoutdf$yD=(extoutdf$ext1-extoutdf$ext2)/NBf
+extoutdf$yT=(extoutdf$ext2)/NBf
+extoutdf$yH=(extoutdf$ext1+extoutdf$ext2)/NBf
+```
+Next conversion and mole fraction is plotted as a function of reaction volume
+
+```{r, tidy=FALSE, fig.width=10, fig.height=5}
+# load library ggplot2 for plotting
+library(ggplot2)
+# load library reshape2 for data reshaping
+library(reshape2)
+# plot of conversion vs volume
+ggplot(extoutdf,aes(x=time,y=conv))+geom_line()+
+  scale_x_continuous(breaks=seq(0,1600,by=200))+xlab('Volume (L)')+ylab('Conversion')+theme_bw(20)
+
+# plot of mole fraction vs volume
+tmp=melt(extoutdf[,c("time","yB","yD","yT","yH")],id.vars=c("time"),variable.name="moleFraction")
+ggplot(tmp,aes(x=time,y=value,color=moleFraction))+geom_line()+
+    scale_x_continuous(breaks=seq(0,1600,by=200))+xlab('Volume (L)')+ylab('moleFraction')+theme_bw(20)
+```
+Info on the session
+```{r, tidy=FALSE}
+print(sessionInfo(),locale=FALSE)
+```   

04_learnR_parfitODE/ABC_data.dat

+   0.000   0.957  -0.031  -0.015
+   0.263   0.557   0.330   0.044
+   0.526   0.342   0.512   0.156
+   0.789   0.224   0.499   0.310
+   1.053   0.123   0.428   0.454
+   1.316   0.079   0.396   0.556
+   1.579   0.035   0.303   0.651
+   1.842   0.029   0.287   0.658
+   2.105   0.025   0.221   0.750
+   2.368   0.017   0.148   0.854
+   2.632  -0.002   0.182   0.845
+   2.895   0.009   0.116   0.893
+   3.158  -0.023   0.079   0.942
+   3.421   0.006   0.078   0.899
+   3.684   0.016   0.059   0.942
+   3.947   0.014   0.036   0.991
+   4.211  -0.009   0.014   0.988
+   4.474  -0.030   0.036   0.941
+   4.737   0.004   0.036   0.971
+   5.000  -0.024   0.028   0.985
Add a comment to this file

04_learnR_parfitODE/ConfRegion_95pct.jpeg

Added
New image
Add a comment to this file

04_learnR_parfitODE/ExpvsPred.jpeg

Added
New image

04_learnR_parfitODE/compAppendix_parest.R

+#
+# Computational Appendix of Book
+# Chemical Reactor Analysis and Design Fundamentals - Rawlings and Ekerdt
+#
+# Example A.5: Estimating rate constants for A->B->C from concentration vs time data
+#
+#
+
+# set working directory
+setwd("~/R/wkspace")
+
+# load libraries
+library(ggplot2) #library for plotting
+library(reshape2) # library for reshaping data (tall-narrow <-> short-wide)
+library(deSolve) # library for solving differential equations
+library(minpack.lm) # library for least squares fit using levenberg-marquart algorithm
+
+#load concentration data
+df=read.table("ABC_data.dat")
+names(df)=c("time","ca","cb","cc")
+
+# plot data
+tmp=melt(df,id.vars=c("time"),variable.name="species",value.name="conc")
+ggplot(data=tmp,aes(x=time,y=conc,color=species))+geom_point(size=3)
+
+# prediction of concentration
+# rate function
+rxnrate=function(t,c,parms){
+  
+  # rate constant passed through a list called parms
+  k1=parms$k1
+  k2=parms$k2
+  
+  # c is the concentration of species
+  
+  # derivatives dc/dt are computed below
+  r=rep(0,length(c))
+  r[1]=-k1*c["A"]  #dcA/dt
+  r[2]=k1*c["A"]-k2*c["B"] #dcB/dt
+  r[3]=k2*c["B"] #dcC/dt
+  
+  # the computed derivatives are returned as a list
+  # order of derivatives needs to be the same as the order of species in c
+  return(list(r))
+  
+}
+
+# predicted concentration for a given parameter set
+cinit=c(A=1,B=0,C=0)
+t=df$time
+parms=list(k1=2,k2=1)
+out=ode(y=cinit,times=t,func=rxnrate,parms=parms)
+head(out)
+
+#plot of predicted concentration
+outdf=data.frame(out)
+tmp=melt(outdf,id.var="time",variable.name="species",value.name="conc")
+ggplot(data=tmp,aes(x=time,y=conc,color=species))+geom_line()
+
+# function that calculates residual sum of squares
+ssq=function(parms){
+ 
+  # inital concentration
+  cinit=c(A=1,B=0,C=0)
+  # time points for which conc is reported
+  # include the points where data is available
+  t=c(seq(0,5,0.1),df$time)
+  t=sort(unique(t))
+  # parameters from the parameter estimation routine
+  k1=parms[1]
+  k2=parms[2]
+  # solve ODE for a given set of parameters
+  out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2))
+  
+  # Filter data that contains time points where data is available
+  outdf=data.frame(out)
+  outdf=outdf[outdf$time %in% df$time,]
+  # Evaluate predicted vs experimental residual
+  preddf=melt(outdf,id.var="time",variable.name="species",value.name="conc")
+  expdf=melt(df,id.var="time",variable.name="species",value.name="conc")
+  ssqres=preddf$conc-expdf$conc
+  
+  # return predicted vs experimental residual
+  return(ssqres)
+  
+}
+
+# parameter fitting using levenberg marquart algorithm
+# initial guess for parameters
+parms=c(k1=0.5,k2=0.5)
+# fitting
+fitval=nls.lm(par=parms,fn=ssq)
+
+# Summary of fit
+summary(fitval)
+# Estimated parameter
+parest=as.list(coef(fitval))
+parest
+# degrees of freedom: # data points - # parameters
+dof=3*nrow(df)-2
+dof
+# mean error
+ms=sqrt(deviance(fitval)/dof)
+ms
+# variance Covariance Matrix
+S=vcov(fitval)
+S
+
+# plot of predicted vs experimental data
+
+# simulated predicted profile at estimated parameter values
+cinit=c(A=1,B=0,C=0)
+t=seq(0,5,0.2)
+parms=as.list(parest)
+out=ode(y=cinit,times=t,func=rxnrate,parms=parms)
+outdf=data.frame(out)
+names(outdf)=c("time","ca_pred","cb_pred","cc_pred")
+
+# Overlay predicted profile with experimental data
+tmppred=melt(outdf,id.var=c("time"),variable.name="species",value.name="conc")
+tmpexp=melt(df,id.var=c("time"),variable.name="species",value.name="conc")
+p=ggplot(data=tmppred,aes(x=time,y=conc,color=species,linetype=species))+geom_line()
+p=p+geom_line(data=tmpexp,aes(x=time,y=conc,color=species,linetype=species))
+p=p+geom_point(data=tmpexp,aes(x=time,y=conc,color=species))
+p=p+scale_linetype_manual(values=c(0,1,0,1,0,1))
+p=p+scale_color_manual(values=rep(c("red","blue","green"),each=2))+theme_bw()
+print(p)
+
+# Get the 95% confidence region
+
+# Inverse of covariance matrix
+Sinv=solve(S)
+
+# draw the confidence region
+# get points for a circle with radius r
+r=sqrt(qf(0.95,2,58)*2)
+theta=seq(0,2*pi,length.out=100)
+z=cbind(r*cos(theta),r*sin(theta))
+# transform points of circle into points of ellipse using 
+# svd of inverse covariance matrix
+Sinv_svd=svd(Sinv)      # inverse of covariance matrix
+xt=t(Sinv_svd$v)%*%diag(1/sqrt(Sinv_svd$d))%*%t(z) # transform from circle to ellispse
+x=t(xt)
+# translate the ellipse so that center is the estimated parameter value
+x=x+matrix(rep(as.numeric(parest),100),nrow=100,byrow=T)
+
+plot(x[,1],x[,2],type="l",xlab="k1",ylab="k2",lwd=2)
+points(parest$k1,parest$k2,pch=20,col="blue",cex=2)
+
+
+# Simulation based estimation of uncertainty
+
+# store original experimental data in a separate dataframe
+dforig=df
+
+# conc profile based on estimated k1 and k2
+cinit=c(A=1,B=0,C=0)
+t=dforig$time
+parms=parest
+out=ode(y=cinit,times=t,func=rxnrate,parms=parms)
+
+outsim=matrix(0,nrow=nrow(dforig),ncol=4)
+outsim[,1]=out[,1]
+
+# number of simulations
+nsim=1000
+
+parsim=matrix(0,nrow=nsim,ncol=2)
+colnames(parsim)=c("k1","k2")
+
+for (i in 1:nsim){
+  
+  # Simulate data set by adding normal random variable with mean 0 and stdev from fit
+  
+  outsim[,2:4]=out[,2:4]+matrix(rnorm(3*nrow(dforig)),nrow=nrow(dforig),ncol=3)*ms
+  df=data.frame(outsim)
+  names(df)=c("time","ca","cb","cc")
+  
+  # get parameter estimate for the simulated dataset
+  parms=as.numeric(parest)
+  fitsim=nls.lm(par=parms,fn=ssq)
+  # store estimated parameters in the ith row
+  parsim[i,]=coef(fitsim)
+  
+  
+}
+
+# plot the parameter estimates from the 1000 simulations
+plot(parsim[,1],parsim[,2],xlab="k1",ylab="k2")
+# overlay the 95% ellipse computed previously
+lines(x[,1],x[,2],col="blue",lwd=2)
+
+# percentage of parameters from simulation within the 95% ellipse
+tmp=rep(0,length.out=nsim)
+for(i in 1:nsim){
+  tmp[i]=(parsim[i,]-as.numeric(parest))%*%Sinv%*%(parsim[i,]-as.numeric(parest))
+}
+sum(tmp <= qf(0.95,2,58)*2)/nsim
+
+# session Info
+sessionInfo()

04_learnR_parfitODE/compAppendix_parest.html

+<!DOCTYPE html>
+<!-- saved from url=(0014)about:internet -->
+<html>
+<head>
+<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
+
+<title>compAppendix_parest.R</title>
+
+<style type="text/css">
+body, td {
+   font-family: sans-serif;
+   background-color: white;
+   font-size: 12px;
+   margin: 8px;
+}
+
+tt, code, pre {
+   font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
+}
+
+h1 { 
+   font-size:2.2em; 
+}
+
+h2 { 
+   font-size:1.8em; 
+}
+
+h3 { 
+   font-size:1.4em; 
+}
+
+h4 { 
+   font-size:1.0em; 
+}
+
+h5 { 
+   font-size:0.9em; 
+}
+
+h6 { 
+   font-size:0.8em; 
+}
+
+a:visited {
+   color: rgb(50%, 0%, 50%);
+}
+
+pre {	
+   margin-top: 0;
+   max-width: 95%;
+   border: 1px solid #ccc;
+   white-space: pre-wrap;
+}
+
+pre code {
+   display: block; padding: 0.5em;
+}
+
+code.r, code.cpp {
+   background-color: #F8F8F8;
+}
+
+table, td, th {
+  border: none;
+}
+
+blockquote {
+   color:#666666;
+   margin:0;
+   padding-left: 1em;
+   border-left: 0.5em #EEE solid;
+}
+
+hr {
+   height: 0px;
+   border-bottom: none;
+   border-top-width: thin;
+   border-top-style: dotted;
+   border-top-color: #999999;
+}
+
+@media print {
+   * { 
+      background: transparent !important; 
+      color: black !important; 
+      filter:none !important; 
+      -ms-filter: none !important; 
+   }
+
+   body { 
+      font-size:12pt; 
+      max-width:100%; 
+   }
+       
+   a, a:visited { 
+      text-decoration: underline; 
+   }
+
+   hr { 
+      visibility: hidden;
+      page-break-before: always;
+   }
+
+   pre, blockquote { 
+      padding-right: 1em; 
+      page-break-inside: avoid; 
+   }
+
+   tr, img { 
+      page-break-inside: avoid; 
+   }
+
+   img { 
+      max-width: 100% !important; 
+   }
+
+   @page :left { 
+      margin: 15mm 20mm 15mm 10mm; 
+   }
+     
+   @page :right { 
+      margin: 15mm 10mm 15mm 20mm; 
+   }
+
+   p, h2, h3 { 
+      orphans: 3; widows: 3; 
+   }
+
+   h2, h3 { 
+      page-break-after: avoid; 
+   }
+}
+
+</style>
+
+<!-- Styles for R syntax highlighter -->
+<style type="text/css">
+   pre .operator,
+   pre .paren {
+     color: rgb(104, 118, 135)
+   }
+
+   pre .literal {
+     color: rgb(88, 72, 246)
+   }
+
+   pre .number {
+     color: rgb(0, 0, 205);
+   }
+
+   pre .comment {
+     color: rgb(76, 136, 107);
+   }
+
+   pre .keyword {
+     color: rgb(0, 0, 255);
+   }
+
+   pre .identifier {
+     color: rgb(0, 0, 0);
+   }
+
+   pre .string {
+     color: rgb(3, 106, 7);
+   }
+</style>
+
+<!-- R syntax highlighter -->
+<script type="text/javascript">
+var hljs=new function(){function m(p){return p.replace(/&/gm,"&amp;").replace(/</gm,"&lt;")}function f(r,q,p){return RegExp(q,"m"+(r.cI?"i":"")+(p?"g":""))}function b(r){for(var p=0;p<r.childNodes.length;p++){var q=r.childNodes[p];if(q.nodeName=="CODE"){return q}if(!(q.nodeType==3&&q.nodeValue.match(/\s+/))){break}}}function h(t,s){var p="";for(var r=0;r<t.childNodes.length;r++){if(t.childNodes[r].nodeType==3){var q=t.childNodes[r].nodeValue;if(s){q=q.replace(/\n/g,"")}p+=q}else{if(t.childNodes[r].nodeName=="BR"){p+="\n"}else{p+=h(t.childNodes[r])}}}if(/MSIE [678]/.test(navigator.userAgent)){p=p.replace(/\r/g,"\n")}return p}function a(s){var r=s.className.split(/\s+/);r=r.concat(s.parentNode.className.split(/\s+/));for(var q=0;q<r.length;q++){var p=r[q].replace(/^language-/,"");if(e[p]){return p}}}function c(q){var p=[];(function(s,t){for(var r=0;r<s.childNodes.length;r++){if(s.childNodes[r].nodeType==3){t+=s.childNodes[r].nodeValue.length}else{if(s.childNodes[r].nodeName=="BR"){t+=1}else{if(s.childNodes[r].nodeType==1){p.push({event:"start",offset:t,node:s.childNodes[r]});t=arguments.callee(s.childNodes[r],t);p.push({event:"stop",offset:t,node:s.childNodes[r]})}}}}return t})(q,0);return p}function k(y,w,x){var q=0;var z="";var s=[];function u(){if(y.length&&w.length){if(y[0].offset!=w[0].offset){return(y[0].offset<w[0].offset)?y:w}else{return w[0].event=="start"?y:w}}else{return y.length?y:w}}function t(D){var A="<"+D.nodeName.toLowerCase();for(var B=0;B<D.attributes.length;B++){var C=D.attributes[B];A+=" "+C.nodeName.toLowerCase();if(C.value!==undefined&&C.value!==false&&C.value!==null){A+='="'+m(C.value)+'"'}}return A+">"}while(y.length||w.length){var v=u().splice(0,1)[0];z+=m(x.substr(q,v.offset-q));q=v.offset;if(v.event=="start"){z+=t(v.node);s.push(v.node)}else{if(v.event=="stop"){var p,r=s.length;do{r--;p=s[r];z+=("</"+p.nodeName.toLowerCase()+">")}while(p!=v.node);s.splice(r,1);while(r<s.length){z+=t(s[r]);r++}}}}return z+m(x.substr(q))}function j(){function q(x,y,v){if(x.compiled){return}var u;var s=[];if(x.k){x.lR=f(y,x.l||hljs.IR,true);for(var w in x.k){if(!x.k.hasOwnProperty(w)){continue}if(x.k[w] instanceof Object){u=x.k[w]}else{u=x.k;w="keyword"}for(var r in u){if(!u.hasOwnProperty(r)){continue}x.k[r]=[w,u[r]];s.push(r)}}}if(!v){if(x.bWK){x.b="\\b("+s.join("|")+")\\s"}x.bR=f(y,x.b?x.b:"\\B|\\b");if(!x.e&&!x.eW){x.e="\\B|\\b"}if(x.e){x.eR=f(y,x.e)}}if(x.i){x.iR=f(y,x.i)}if(x.r===undefined){x.r=1}if(!x.c){x.c=[]}x.compiled=true;for(var t=0;t<x.c.length;t++){if(x.c[t]=="self"){x.c[t]=x}q(x.c[t],y,false)}if(x.starts){q(x.starts,y,false)}}for(var p in e){if(!e.hasOwnProperty(p)){continue}q(e[p].dM,e[p],true)}}function d(B,C){if(!j.called){j();j.called=true}function q(r,M){for(var L=0;L<M.c.length;L++){if((M.c[L].bR.exec(r)||[null])[0]==r){return M.c[L]}}}function v(L,r){if(D[L].e&&D[L].eR.test(r)){return 1}if(D[L].eW){var M=v(L-1,r);return M?M+1:0}return 0}function w(r,L){return L.i&&L.iR.test(r)}function K(N,O){var M=[];for(var L=0;L<N.c.length;L++){M.push(N.c[L].b)}var r=D.length-1;do{if(D[r].e){M.push(D[r].e)}r--}while(D[r+1].eW);if(N.i){M.push(N.i)}return f(O,M.join("|"),true)}function p(M,L){var N=D[D.length-1];if(!N.t){N.t=K(N,E)}N.t.lastIndex=L;var r=N.t.exec(M);return r?[M.substr(L,r.index-L),r[0],false]:[M.substr(L),"",true]}function z(N,r){var L=E.cI?r[0].toLowerCase():r[0];var M=N.k[L];if(M&&M instanceof Array){return M}return false}function F(L,P){L=m(L);if(!P.k){return L}var r="";var O=0;P.lR.lastIndex=0;var M=P.lR.exec(L);while(M){r+=L.substr(O,M.index-O);var N=z(P,M);if(N){x+=N[1];r+='<span class="'+N[0]+'">'+M[0]+"</span>"}else{r+=M[0]}O=P.lR.lastIndex;M=P.lR.exec(L)}return r+L.substr(O,L.length-O)}function J(L,M){if(M.sL&&e[M.sL]){var r=d(M.sL,L);x+=r.keyword_count;return r.value}else{return F(L,M)}}function I(M,r){var L=M.cN?'<span class="'+M.cN+'">':"";if(M.rB){y+=L;M.buffer=""}else{if(M.eB){y+=m(r)+L;M.buffer=""}else{y+=L;M.buffer=r}}D.push(M);A+=M.r}function G(N,M,Q){var R=D[D.length-1];if(Q){y+=J(R.buffer+N,R);return false}var P=q(M,R);if(P){y+=J(R.buffer+N,R);I(P,M);return P.rB}var L=v(D.length-1,M);if(L){var O=R.cN?"</span>":"";if(R.rE){y+=J(R.buffer+N,R)+O}else{if(R.eE){y+=J(R.buffer+N,R)+O+m(M)}else{y+=J(R.buffer+N+M,R)+O}}while(L>1){O=D[D.length-2].cN?"</span>":"";y+=O;L--;D.length--}var r=D[D.length-1];D.length--;D[D.length-1].buffer="";if(r.starts){I(r.starts,"")}return R.rE}if(w(M,R)){throw"Illegal"}}var E=e[B];var D=[E.dM];var A=0;var x=0;var y="";try{var s,u=0;E.dM.buffer="";do{s=p(C,u);var t=G(s[0],s[1],s[2]);u+=s[0].length;if(!t){u+=s[1].length}}while(!s[2]);if(D.length>1){throw"Illegal"}return{r:A,keyword_count:x,value:y}}catch(H){if(H=="Illegal"){return{r:0,keyword_count:0,value:m(C)}}else{throw H}}}function g(t){var p={keyword_count:0,r:0,value:m(t)};var r=p;for(var q in e){if(!e.hasOwnProperty(q)){continue}var s=d(q,t);s.language=q;if(s.keyword_count+s.r>r.keyword_count+r.r){r=s}if(s.keyword_count+s.r>p.keyword_count+p.r){r=p;p=s}}if(r.language){p.second_best=r}return p}function i(r,q,p){if(q){r=r.replace(/^((<[^>]+>|\t)+)/gm,function(t,w,v,u){return w.replace(/\t/g,q)})}if(p){r=r.replace(/\n/g,"<br>")}return r}function n(t,w,r){var x=h(t,r);var v=a(t);var y,s;if(v){y=d(v,x)}else{return}var q=c(t);if(q.length){s=document.createElement("pre");s.innerHTML=y.value;y.value=k(q,c(s),x)}y.value=i(y.value,w,r);var u=t.className;if(!u.match("(\\s|^)(language-)?"+v+"(\\s|$)")){u=u?(u+" "+v):v}if(/MSIE [678]/.test(navigator.userAgent)&&t.tagName=="CODE"&&t.parentNode.tagName=="PRE"){s=t.parentNode;var p=document.createElement("div");p.innerHTML="<pre><code>"+y.value+"</code></pre>";t=p.firstChild.firstChild;p.firstChild.cN=s.cN;s.parentNode.replaceChild(p.firstChild,s)}else{t.innerHTML=y.value}t.className=u;t.result={language:v,kw:y.keyword_count,re:y.r};if(y.second_best){t.second_best={language:y.second_best.language,kw:y.second_best.keyword_count,re:y.second_best.r}}}function o(){if(o.called){return}o.called=true;var r=document.getElementsByTagName("pre");for(var p=0;p<r.length;p++){var q=b(r[p]);if(q){n(q,hljs.tabReplace)}}}function l(){if(window.addEventListener){window.addEventListener("DOMContentLoaded",o,false);window.addEventListener("load",o,false)}else{if(window.attachEvent){window.attachEvent("onload",o)}else{window.onload=o}}}var e={};this.LANGUAGES=e;this.highlight=d;this.highlightAuto=g;this.fixMarkup=i;this.highlightBlock=n;this.initHighlighting=o;this.initHighlightingOnLoad=l;this.IR="[a-zA-Z][a-zA-Z0-9_]*";this.UIR="[a-zA-Z_][a-zA-Z0-9_]*";this.NR="\\b\\d+(\\.\\d+)?";this.CNR="\\b(0[xX][a-fA-F0-9]+|(\\d+(\\.\\d*)?|\\.\\d+)([eE][-+]?\\d+)?)";this.BNR="\\b(0b[01]+)";this.RSR="!|!=|!==|%|%=|&|&&|&=|\\*|\\*=|\\+|\\+=|,|\\.|-|-=|/|/=|:|;|<|<<|<<=|<=|=|==|===|>|>=|>>|>>=|>>>|>>>=|\\?|\\[|\\{|\\(|\\^|\\^=|\\||\\|=|\\|\\||~";this.ER="(?![\\s\\S])";this.BE={b:"\\\\.",r:0};this.ASM={cN:"string",b:"'",e:"'",i:"\\n",c:[this.BE],r:0};this.QSM={cN:"string",b:'"',e:'"',i:"\\n",c:[this.BE],r:0};this.CLCM={cN:"comment",b:"//",e:"$"};this.CBLCLM={cN:"comment",b:"/\\*",e:"\\*/"};this.HCM={cN:"comment",b:"#",e:"$"};this.NM={cN:"number",b:this.NR,r:0};this.CNM={cN:"number",b:this.CNR,r:0};this.BNM={cN:"number",b:this.BNR,r:0};this.inherit=function(r,s){var p={};for(var q in r){p[q]=r[q]}if(s){for(var q in s){p[q]=s[q]}}return p}}();hljs.LANGUAGES.cpp=function(){var a={keyword:{"false":1,"int":1,"float":1,"while":1,"private":1,"char":1,"catch":1,"export":1,virtual:1,operator:2,sizeof:2,dynamic_cast:2,typedef:2,const_cast:2,"const":1,struct:1,"for":1,static_cast:2,union:1,namespace:1,unsigned:1,"long":1,"throw":1,"volatile":2,"static":1,"protected":1,bool:1,template:1,mutable:1,"if":1,"public":1,friend:2,"do":1,"return":1,"goto":1,auto:1,"void":2,"enum":1,"else":1,"break":1,"new":1,extern:1,using:1,"true":1,"class":1,asm:1,"case":1,typeid:1,"short":1,reinterpret_cast:2,"default":1,"double":1,register:1,explicit:1,signed:1,typename:1,"try":1,"this":1,"switch":1,"continue":1,wchar_t:1,inline:1,"delete":1,alignof:1,char16_t:1,char32_t:1,constexpr:1,decltype:1,noexcept:1,nullptr:1,static_assert:1,thread_local:1,restrict:1,_Bool:1,complex:1},built_in:{std:1,string:1,cin:1,cout:1,cerr:1,clog:1,stringstream:1,istringstream:1,ostringstream:1,auto_ptr:1,deque:1,list:1,queue:1,stack:1,vector:1,map:1,set:1,bitset:1,multiset:1,multimap:1,unordered_set:1,unordered_map:1,unordered_multiset:1,unordered_multimap:1,array:1,shared_ptr:1}};return{dM:{k:a,i:"</",c:[hljs.CLCM,hljs.CBLCLM,hljs.QSM,{cN:"string",b:"'\\\\?.",e:"'",i:"."},{cN:"number",b:"\\b(\\d+(\\.\\d*)?|\\.\\d+)(u|U|l|L|ul|UL|f|F)"},hljs.CNM,{cN:"preprocessor",b:"#",e:"$"},{cN:"stl_container",b:"\\b(deque|list|queue|stack|vector|map|set|bitset|multiset|multimap|unordered_map|unordered_set|unordered_multiset|unordered_multimap|array)\\s*<",e:">",k:a,r:10,c:["self"]}]}}}();hljs.LANGUAGES.r={dM:{c:[hljs.HCM,{cN:"number",b:"\\b0[xX][0-9a-fA-F]+[Li]?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+(?:[eE][+\\-]?\\d*)?L\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+\\.(?!\\d)(?:i\\b)?",e:hljs.IMMEDIATE_RE,r:1},{cN:"number",b:"\\b\\d+(?:\\.\\d*)?(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\.\\d+(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"keyword",b:"(?:tryCatch|library|setGeneric|setGroupGeneric)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\.",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\d+(?![\\w.])",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\b(?:function)",e:hljs.IMMEDIATE_RE,r:2},{cN:"keyword",b:"(?:if|in|break|next|repeat|else|for|return|switch|while|try|stop|warning|require|attach|detach|source|setMethod|setClass)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"literal",b:"(?:NA|NA_integer_|NA_real_|NA_character_|NA_complex_)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"literal",b:"(?:NULL|TRUE|FALSE|T|F|Inf|NaN)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"identifier",b:"[a-zA-Z.][a-zA-Z0-9._]*\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"<\\-(?!\\s*\\d)",e:hljs.IMMEDIATE_RE,r:2},{cN:"operator",b:"\\->|<\\-",e:hljs.IMMEDIATE_RE,r:1},{cN:"operator",b:"%%|~",e:hljs.IMMEDIATE_RE},{cN:"operator",b:">=|<=|==|!=|\\|\\||&&|=|\\+|\\-|\\*|/|\\^|>|<|!|&|\\||\\$|:",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"%",e:"%",i:"\\n",r:1},{cN:"identifier",b:"`",e:"`",r:0},{cN:"string",b:'"',e:'"',c:[hljs.BE],r:0},{cN:"string",b:"'",e:"'",c:[hljs.BE],r:0},{cN:"paren",b:"[[({\\])}]",e:hljs.IMMEDIATE_RE,r:0}]}};
+hljs.initHighlightingOnLoad();
+</script>
+
+
+
+
+</head>
+
+<body>
+<!-- Automatically generated by RStudio [12861c30b10411e1afa60800200c9a66] -->
+
+<h3>compAppendix_parest.R</h3>
+
+<p>admin &mdash; <em>Jun 30, 2013, 10:43 AM</em></p>
+
+<pre><code class="r">#
+# Computational Appendix of Book
+# Chemical Reactor Analysis and Design Fundamentals - Rawlings and Ekerdt
+#
+# Example A.5: Estimating rate constants for A-&gt;B-&gt;C from concentration vs time data
+#
+#
+
+# set working directory
+setwd(&quot;~/R/wkspace&quot;)
+
+# load libraries
+library(ggplot2) #library for plotting
+library(reshape2) # library for reshaping data (tall-narrow &lt;-&gt; short-wide)
+library(deSolve) # library for solving differential equations
+library(minpack.lm) # library for least squares fit using levenberg-marquart algorithm
+
+#load concentration data
+df=read.table(&quot;ABC_data.dat&quot;)
+names(df)=c(&quot;time&quot;,&quot;ca&quot;,&quot;cb&quot;,&quot;cc&quot;)
+
+# plot data
+tmp=melt(df,id.vars=c(&quot;time&quot;),variable.name=&quot;species&quot;,value.name=&quot;conc&quot;)
+ggplot(data=tmp,aes(x=time,y=conc,color=species))+geom_point(size=3)
+</code></pre>
+
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAA21BMVEUAAAAAADoAAGYAOjoAOpAAZrYAujg6AAA6ADo6AGY6OgA6Ojo6OpA6ZmY6ZrY6kNthnP9mAABmADpmAGZmZgBmZjpmkNtmtv9/f39/f5V/f6t/lcF/q9aQOgCQOjqQkGaQ2/+Vf3+VlcGVq9aVweurf6urlZWrlcGr1v+2ZgC225C22/+2///BlX/BlZXBlavBq8HBwdbB6//Wq3/Wq5XW///bkDrb25Db///l5eXrwZXr1qvr///y8vL4dm3/tmb/1qv/25D/27b/68H//7b//9b//9v//+v///8fExWrAAAACXBIWXMAAAsSAAALEgHS3X78AAAS+0lEQVR4nO3dC3/T1h3GccMgHds6oN1Wc+tgCxA6LmvSZJQFx6wk1vt/RZMs33Txkf7SkfTI53c+FGIRnurJ1zo+thVlEjGCHJOhd4AxzAA+0AF8oAP4QAfwgQ7gAx3ABzqAD3QAH+gAPtABfKAD+EAH8IGORvC/5kdxS358rvoEDxEqGZ3sBvD6GcBbqoqgyVYBXj8DeEtVETTZKsDrZwBvqSqCJlsFeP0M4C1VRdBkqwCvnwG8paoImmwV4PUzgLdUFUGTrQK8fgbwlqoiaLJVeoS/fnm+/PNs+uh0/Tvwhw8/nz5ews+fXG5+AV8nY9zwi58XH5bwn46TYz/9Pb45jYfvnTig8eDB0HtQczim+jX8SfJR+vvqb7q4j4scrS0zHsRj5Ef8Fj57xAPvyHiwHFFum5fd6B+ex/j6GSXwyznAvBvFf9Mz/PXrS1b1hozCVJ/eFay7UfJveoTfPwo7Cvx6JGCt4cv+EfD6Ga2neuBtESoZ7Rd3TPWmCJUMc0TxvjHo4g74hhnWiLJHA96kMUSoZBgjStd/wBsiVDKA31e1bIigMdUD3zjDw+JuJPCzmbVqcYigyd6HFeFns2p54FvuhiD8bFZDHviWuwG8fkYo8Ez1feyGIjyLux52QxJe9qvVJqPszZV8RvnnVOxG/I+A9xrhNaP07dRcxp7Pce9G8o+A9xrhM6P8BIpsxr7Pce6GI9iZAXw/GcAHCs9UHyo8i7tQ4Ztn8HTOUlUETbYK8PoZwFuqiqDJVgFePwN4S1URNNkqwOtnAG+pKoImWwV4/QzgLVVF0GSrAK+fAbylqgiaK2P1yjzwXqvqw6/fiwPea1V5+M0b6cB7rTo0/PodVuDDgt+cU8FUHxT89oQoFnfAN90P4C1V9af62vsBvKWq/OKu/n4Ab6k6NLzHDOAtVUXQZKsAr58BvKWqCJpsFQn4z/lR3GIeHiJUMjrZDQn4Lu7jIkerbBXg9TOAt1QVQZOtArx+BvCWql2ibb+pFfiQ4He+jR34gOB3r0EBPPDd7Qfwlqr9TvX7r0KjWgX4BhmFxZ3j+kOqVYD3kOG68pRqFeA9ZAAfKDxTfajwLO5Che80A3hLVRE02SrA62cAb6kqgiZbBXj9DOAtVUXQZKsAr58BvKWqCJpsFeD1M4C3VBVBk60CvH4G8JaqImiyVYDXzwDeUlUETbYK8PoZwFuqiqDJVgFePwN4S1URNNkqwOtnAG+pKoImWwV4/QzgLVVF0GSrAK+fAbylqgiabBXg9TOAt1QVQZOtArx+BvCWqiJoslWA188A3lJVBE22CvD6GcBbqoqgyVYBvjJj/0UPgK8DfzZ9dJr8+Wkaj+P45uq2OrzjMifA14CfP7mMf6UfL/59uXh/uv27Lqr6+mq5LmwEfA34T8fR9cvz9ONfTqPrVy+mT5OPk+Pf9074HCn80HuhP/bDn0SLDyn8PBafPz5PNqWji/s4U707oz/47RH/y2qWnx+PAZ7FXTv47WP89evLFH0cR/zwGeOGT1f1Cfpv/0xvTtcHPPAD7EZ/8I7RRVURNNkqwOtnAG+pKoImWwV4/QzgLVVF0GSrAK+fAbylqgiabBXg9TOAt1QVQZOtArx+BvCWqiJoslWA188A3lJVBE22CvD6GcBbqoqgyVYBXj8DeEtVETTZKsDrZwBvqSqCJlsFeP0M4C1VRdBkqwCvnwG8paoImmwV4PUzgLdUFUGTrQK8fgbwlqoiaLJVgNfPAN5SVQRNtgrwubH73fDAhwOfuf4F8IcLn7vcRfZSR8AfLHz+AjfAhwFfvJYZU32g8CzugoB3X8sM+MOFd17LDPgDhh9DBvCWqiJoslWA188A3lJVBE22CvD6GcBbqoqgyVYBXj8DeEtVETTZKhLwn/OjuMU8PESoZHSyGxLwXdzHRY5W2SrA62cAb6kqgiZbBXj9DOAtVUXQZKsAr58BvKWqCJpsFeD1M4C3VBVBk60CvH4G8JaqImiyVYDXzwDeUlUETbYK8PoZwFuqiqDJVgFePwN4S1URNNkqwOtnAG+pKoImWwV4/QzgLVVF0GSrAK+fAbylqgiabJWRwM9m1qoiaMC3gp/N8vLNvlqZ618ALw8/mxXkG321sle8AT4U+Nw1roCXh/c01QO/u6HW+PL754PCe1rcMdXvbPA8tJ/OsbjbbsiMm2eTyd3o6tYfJnc+Rm8mk/tRdJFsSY745c30E0YLb4tQyegBPgb+379i+OdfH96/uP3268Mflh//EG9Pb6afAHyvGX1M9fFhffvt1e230Zu78YfxMX4RH/rLO0R6M/0E4HvN6OcxPjnKJz8kR/ySPD7ib57dT4745c30E4DvNaOnx/g7H69u/yX+Pf04eYy/8zGGT2+uNgLfZ0Zvq/or52TuHsD7zwDeUlUETbZKY2Hge8sA3lJVBE22So6g9dXQgPefAbylqgiabJWxwh8dWZsWBvAjhD86qpIXQRsZ/GymDX90VCkvgjYu+OS8B+BlMnqDT8908gJ/NfkhurhVdRaHvSpTfdvdqAG/fLX+TcX78OXwX79LzL/cq3glsEHVontuiwjaqOCzU/2XP32MLr79c3TzY53Tr7LwN8+Sd/Kulm/rnE0fnS43nk2TjzY3m8EXqubnABG0ccFnFndX95O/iA/7yhm7CB99fThZvYE/f3IZ/4o/WLw/3b3pCb7wqC+CNjL43adzyRF/lZx11eCI3xmfjqPrl+fxB9evXkyfbm9O41Ejt2qk8B6CAh1lz+Pjg/32T99Mbv3VdQLGeuyHP4kWHxLp+ePz+MbmZjJ83MeZ6m0ZNeBNo7Cqn6zm+s0hnoz5ceaml6os7kwZ3cLvnKa1eVCfHyeHv+/H+C4iVDLGCP/99oncchl//foyWdUfR95X9R1EqGSMED56c7/OP+qiqgiabBWzrAl++Wyu4oRs4AfZjRyB7yO+3uiiqgiabBXg9TPGCJ+ch+8+ER/4Q4S/eZYs7i6q5LuoKoImW6Vb+PTp3O6TOuDtGT3CHx1xxOtk9AefvO6dhf/y50bwPMbLVimDT9/p8gJfb3RRVQRNtkoN+OQMnP98+7DJGTjDwmcufNIgIjj47FSfnoFz7230psHbsl2dc1fnq5W91FGDiPDgM4u75Rk4X/62PhXHBN/dOXfVX63cxc0aRAQIv/t0bnkGzu/i3xoc8bvn3AHfNGOo5/HLM3D++E2jx/jtOXe9wzPVOzNqwJsGizv/GcBbqoqgyVYBXj9jhCdiAC9bJUfAES+YAbylqgiabJVDhq/4hlrgDxS+6lupgT9M+MqLJwBfhH/wAHg/uzEu+OTFztHDM9U7M8rg07c3Rg/P4s6VUQN+eSmUixpnUKnBt41QyRhoqk9PxEh+s78fD7xslVL4zOJuef5FrZMwgO8mY+gTMa7qvCEPvP+MIU/EeHtR43wK4LvJ4AUcS1URNNkqwOtnAG+pWhqReWYPfDDw2dfygPc5lOFzr94DD3z9DOtuAC8Az1S/uyEkeBZ3OxuCgrdFqGQAb6kqgiZbBXj9DOAtVUXQZKsAr59xuPD51ws/F7eYh4cIlYxOdkMCvov7uMjRKlsFeP0M4C1VRdBkqwCvnwG8paoImmwV4PUzgLdUFUGTrQK8fgbwlqoiaLJVxgo/m1mbFgbwI4RPftC5sWlhAD8++NmsUl4EDXjgG2cAvx1M9W13YyTwhYtbsLhruRvjgC9ezoancy13YxTwJRewAr7lbgCvnxEuPFO9/90YB3xxcQd8y90YCby9qgiabBXg9TOAd1bN/eghETTgu4bP/7AxETTgO4Yv/HhBETTgReGdFzoGXh++4VTvvrQ58COAb7S4q/hhBsCPAb5JBPCBwjPVhwrP4i5U+MLYvScAHw58Zu4HPhj47GoPeOCb7gbw44Bnqg8VnsVdqPCCGcBbqoqgyVYBXj8DeEtVETTZKsDrZwBvqSqCJlsFeP0M4C1VRdBkqwCvnwG8paoImmyV/uDPpo9Okz8X76bTJ5fxzdVt4AfZjd7g508u41/JB0/jO8HJ4v3p9u+6qCqCJlulN/hPx9H1y/P1jZPrVy+mT5MPp/HwvROM/sd++JNo8WEFHx/088fnyaZ0dHEfFzlaZav0B7894s+ervyPga+TMW74zWP84t3JCp0jvl7GuOHTVf3168uz5FH9OFnVrw944AfYjf7gHaOLqiJoslWA188A3lk1d61LETTgu4bPX91WBA34juEL17P299VyfXsd8IcL7/yGWuCHhu9sqnd/Cz3wg8N3tbgDXh3efwRTfdDwLO5ChR88A3hLVRE02SrA62cAb6kqgiZbBXj9DOAtVUXQZKsAr58BvKWqCJpsFeD1M4C3VO0Xbf+re6pVgPeR4Xg9X7UK8B4yXO/gqVYB3kMG8IHCM9WHCs/iLlT4TjOAt1QVQZOtArx+BvCWqiJoslWA188AfmfkfrxgjaaFAfwI4fM/ULRG08IAfnzwhR8hXKNpYQAPvDFj+3IN8P3CDzvV77xAC3zP8EMu7nbfkgG+b/gBn84BHyg8U32o8CzuQoXvOQN4S1URNNkqwOtnAG+pKoImW+WQ4XNXxTFHAD9O+Px1sMwRwHcN/zk/ilvMI0qvfNcuo/1ueKnSPqKYIQHfxX08Klzy0BzRajfWz/U54vuGH3aq37y6B3zv8EMu7rav5wPfP3zbCOCBN2cw1QcKz+IuVHiPGcBbqoqgyVYBXj8DeEtVETTZKsDrZwBvqSqCJlsFeP0M4C1VRdBkqwCvnwG8pWppROZtm57hi9dHWm0B3mvVsojsG7X9wheviLbeArzXqiURuVMzeoUvXgNxswV4r1WBd2eEBM9Uv7shJHgWdzsbgoK3RahkAG+pKoImWwV4pYzyayEDb6k6Rvg9Vz8H3lJ1hPD7ft4B8JaqwLszgBfKYKoPFJ7FXWFL4bJ2hwnf426MA754IUvgW+7GKOBLLl0LfMvdCB3e+R21AvDJo35FRPwpnwtbKndjFPDdTfXu76EfHn65zndHJJ/yubClcjfGAd/V4q7iqhmDw6fP7J0Rjjf6nbsxEvjC6Ax+9+bhwJdMAUHDF6f6zIZotanFbmhM9WVTQNjwedXsFBCtNzXfDYnFXdkcEDp8bhTh3cuAcTydA746ojDVl8Fvb48Dnqm+RkRhcVfmvtniGb58avHw1WBxZ4nYs7jbnQP8wpc8qCQbshFlb+Qc6it3JWPAl2w7gy99chlvyESUvnULvNeIfRm1pvr1p7SCT7fsRpSfrAG814i9GTUWdxtCx36sPsMx1QNfGOLvzm0P3v0Z689wLe6aTfUla0TgDRHdwm8+xbkfTRZ3ZS88DAh/Nn10uvPB5uZBwJcfqs6MevC23cgGOzN6g58/uYx/bT7Y3DwI+D3Pw9wZJVN9y93Y5grBfzqOrl+ebz7Y3JzGw/dO9D1WazD7P+tgX6J0d7pJdoz98CfR4sP55oPNzWS0vY+XjT6PeNfr+d6rlP2PlBd3+474Q4B3vYPnu0rp/0r56ZzlMb7yZ4SrwTves/cMXz65KMOny/jr15fVq/riGXbFIQbfacbI4R0ju0sl59QWR6jw45vqg4ff+2jgf3FXmfGrJPxhTvX71399VCl5KU8QfoSLu8oMxzO+HqqUvaKjCD+6p3PVGR3Cx7EVEWX/c0n48u8WzoyxwXc31SfBBwK/5/oAmTE6eG+Lu9xwvWqY+aR8qB78viuCZMb44JtntIcfx+IOeFNEjam+LEMQ/kCn+uYZrRd3ZRmK8Ae5uGuR0cluSMLLfrWGyQDeUlUETbYK8PoZwFuqiqDJVgFePwN4S1URNNkqwOtnAG+pKoImWwV4/QzgLVVF0GSrAK+fAbylqgiabBXg9TOAt1QVQZOtArx+BvCWqiJoslWA1884XPjC8HCpBB9XWxDJENkN9wDef4bIbrgH8P4zRHbDPYD3nyGyG+7hB54xugF8oAP4QAfwgQ7gAx0+4HcviNV0bC+i13Qs3k2nmwuyNR1n0+lJy4jllQHb70b7L6l7eIDPXAKvacb0cVv4+dP4C3bSOqP9PXA+bQu/eN+xeuQFPnPRy2Zj8fP2eqmtduWkdURr+N/+8d+28NevXkyftsyoGj7gdy9z23R4gZ+3/2qdtT1c4zvOvC38PJ7+PNyHnUPjiPcDf+blKGn5CP0pubx36wf5GN9DhmuIPMb7mDPenbRMSL/a7Zdm7Y/4Yy+PWs6hsqpvD3/m41CLQ9pPG+0P1jMvk4Zz8Dw+0AF8oAP4QAfwgQ7gAx0BwH+59zb+NfReqI0w4IfeBcFx+PBfH05u/3Tv7c2Pf59M7l/F/6WbQr8zHD78aqq/eXY3+vLN3eXx/+Z+dHF36N0aeIQD/+PzKPnv63fP41/R1+8DP+SDhH84mUxuPR96v4YdQcKHfrQnI0T45DH+6s7Hofdr2BEA/M2zdFW/hY/n+tBn+hDgGWUD+EAH8IEO4AMdwAc6gA90AB/oAD7QAXygA/hAx/8Bq3qc0q9HjawAAAAASUVORK5CYII=" alt="plot of chunk unnamed-chunk-1"/> </p>
+
+<pre><code class="r">
+# prediction of concentration
+# rate function
+rxnrate=function(t,c,parms){
+
+  # rate constant passed through a list called parms
+  k1=parms$k1
+  k2=parms$k2
+
+  # c is the concentration of species
+
+  # derivatives dc/dt are computed below
+  r=rep(0,length(c))
+  r[1]=-k1*c[&quot;A&quot;]  #dcA/dt
+  r[2]=k1*c[&quot;A&quot;]-k2*c[&quot;B&quot;] #dcB/dt
+  r[3]=k2*c[&quot;B&quot;] #dcC/dt
+
+  # the computed derivatives are returned as a list
+  # order of derivatives needs to be the same as the order of species in c
+  return(list(r))
+
+}
+
+# predicted concentration for a given parameter set
+cinit=c(A=1,B=0,C=0)
+t=df$time
+parms=list(k1=2,k2=1)
+out=ode(y=cinit,times=t,func=rxnrate,parms=parms)
+head(out)
+</code></pre>
+
+<pre><code>      time       A      B       C
+[1,] 0.000 1.00000 0.0000 0.00000
+[2,] 0.263 0.59096 0.3556 0.05348
+[3,] 0.526 0.34924 0.4834 0.16731
+[4,] 0.789 0.20639 0.4958 0.29779
+[5,] 1.053 0.12172 0.4543 0.42395
+[6,] 1.316 0.07193 0.3925 0.53552
+</code></pre>
+
+<pre><code class="r">
+#plot of predicted concentration
+outdf=data.frame(out)
+tmp=melt(outdf,id.var=&quot;time&quot;,variable.name=&quot;species&quot;,value.name=&quot;conc&quot;)
+ggplot(data=tmp,aes(x=time,y=conc,color=species))+geom_line()
+</code></pre>
+
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAA0lBMVEUAAAAAADoAAGYAOpAAZrYAujg6AAA6ADo6AGY6OgA6OpA6kNthnP9mAABmADpmAGZmOpBmZgBmZmZmtv9/f39/f5V/f6t/lcF/q9aQOgCQOjqQZpCQkGaQ2/+Vf3+VlcGVq9aVweurf6urlZWrlcGr1v+2ZgC225C22/+2///BlX/BlZXBlavBq8HBwdbB6//Wq3/Wq5XW///bkDrb/7bb///l5eXrwZXr1qvr///y8vL4dm3/tmb/1qv/25D/27b/68H//7b//9b//9v//+v///9BdM7AAAAACXBIWXMAAAsSAAALEgHS3X78AAAWN0lEQVR4nO3dC3vTyBUGYC+g7baIblvMbVu2DRBaLt1kQ3GaOKabxP7/f6kay5Y0mrvmorl85wlgO+b4zLweWZZG0mKHKDIWcxeAmCcAX2gAvtAAfKEB+EID8IUG4AsNwBcagC80AF9oAL7QAHyhAfhCYxL8ahzsI+O4Vj3BQYpYcngpA/Dx5wC8SVMjQYu2KYCPPwfgTZoaCVq0TQF8/DkAb9LUSNCibQrg488BeJOmRoIWbVMAH38OwJs0NRK0aJsC+PhzAN6kqZGgRdsUwMefA/AmTY0ELdqmAD7+HIA3aWokaNE2BfDx50gd/u7Nxf7f8+Wzs+PfgM8ffrN8voffvLjqfgCvkyNt+O0v2897+MsTMvbbv5u7yyaYJ69dl4XwHZJF/RH+lNxq/z78hnl/ru3f45GMVoz4Hp4e8YAvBV7nMx7w3ssIDH/37kpnrR7w3ssICC8OtkylPOAtywB8/DlsU9RNAN4gRSw5Jqaoj8HNAfj4c5iloLzFOQAffw6tFHxvcQ7Ax59DlkLhLc4B+Phz8FLQ3lirl7TUPEUsOagU3AEOeElLzVPEkqNLIVykA17SUvMUseQgKeSf4/nAKzfalgSvXncDvKSl5iliyMHd6uaiDMDHm0O41c1FGYCPM8dw8Q54k6amDD/6TAe8SVNTheesyQHepKlJwvNX3wFv0tT04PV3sLgoA/Bx5JB9VS8KXiWfFbxiAw3gTZqaDPyEXaouygD8rDmIOiZbUmUWAH8Y64CnyswevlvCA54qM2/44Qc74KkyM4Yfrc4BniozW3hmJR7wVJl5wvO+uwGeLlMunyQ8/xs74Okys4MXbagBPF1mZvDi7XOAp8vMCl62WRbwdJk5wUu3xgOeLjMfeMVeGMDTZeYCr9z5Bni6zDzgleyAH5eZBXztog7AmzQ1Bngy3AFvCi+XTwDe2WEwgDdp6uzwtbM6AG/S1Jnhu5U6wBcF36/LA74k+MF3OMCXA099dwd8MfCKMwua1wF4k6bOBT/eVAf4MuCZLbSAl8T1ONpH1szjBsEmDZGj9lGHl6ZEAS94f8qGfJQjnrdHBiM+f3jujjjAZw/P3wEL+MzhPZ7OAvAmTQ0ML5xvAfis4b2exwTwJk0NCi+ZXwX4jOFl8+oAbw4vk48I3vuZiwBv0tRg8IpptIDPFD7AKasAb9LUQPDKafOAzxJefbgE4HOEr4McjQN4k6aGgNc4OSHgM4SvwxyGBXijpnqH3399B3xx8LWDHJp1lAcvkZ8bfvpFXc3rALxJU/3Ca5+HFvBZwXdf3wFfFLzBCYgBnxH8YHMd4AuCH26mBTzg/dUBeJOmeoM3O9c84HOBp/fHAd4PvFh+LvjRfljAFwJvfHUJwGcBz8y7AHwR8BMuKwJ4wE+uA/AmTfXQW5wJdoAvAJ43sRLw+cNPu4IU4FOH50+kBrwneKF86N6aeukwwKcNLzpwAvCA91cH4E2a6rS3pl8sEPApw4uPkAN8zvA2V4kEfLrwsiNiAQ94f3UA3qSpznrL7rqwgE8V3vK6sICfCC+SB7xlGYCX51Cc6wTwavjz5bMz8u/lsomT5u7hftTw1leCBvzmxVXz097e/vtq++ms/52szFnhlec2ArwS/vJkd/fmor3969nu7u3Py5fkNhn/soRrl9WZRj3ni6cVYvjT3fZzC79pxDfPL8hDbcjen3OOeAfXfseIH4z4Xw9L+c1J3PC1lx43z5E2fP8Zf/fuqkWPfcQD3gV8u1ZP0H/7R3t3eRzwkcLXfnrcPEfi8JKQlsmX999btaceN88BeJOmAl6eA/CCqB3kcFGHrzIAz4/aQQ4XdXgrA/D8AHyZ8LWDHC7q8FcG4HlRO8jhog6PZQCeE8dNtYAHvHEOF3X4LCN+eL68z97q9s0Avij4fp8c4AFvmMNFHX7LAPw4BjvhAV8Q/HDyBeABb5TDRR2+ywA8HdRsK8AD3iSHizq8lwF4KujplYAHvEEOF3X4LyMBeK68n94azacGPOD1c7ioI0AZgB/E+AAKwANeO4eLOkKUAfg+mCOmAA943Rwu6ghSBuC7YA+RBDzgNXO4qCNMGYA/BueYaMDPC8+Td95bvGPhAQ94rRwu6ghVBuDb4J78AvCA18nhoo5gZQB+H/yz3QAe8Bo5XNQRrgzAkxCc3grwgFfncFFHwDIAvxKfzw7wM8Nz5AFvWQbgJSewBDzgVTlc1BG0DMBLzlgLeMArcrioI2wZgJecohrwgJfncFFH4DKKh5edkx7wgJfmcFFH6DKigL8ex/iRNfMMZbBJuVE7yOGiDs8p2BxRwKvfn8yQdzVMpFcfwYjPFl5+1RnAA16Sw0Ud4csoG15xmSnAA16cw0UdM5RRNLzqunKAB7wwh4s65iijZHjlhSQBD3hRDuUzAG8Fz8g76C31lWMBD3hBAD49+DoIGuABPzkH4E2aquqtOgwa4AE/OQfgTZqq6K06EBrgLeHH8oC3LAPw03PolAH4uOBrnRSAB7woAJ8WfK2VAvC7b797DXhe5A6vH97gR/JWvVXrpSgI/v7VYvFod/Pd7xcPv+7eLxaPd7sv5BEy4vd32ycAXruMROAb4P/9q4F/ffvk8ZcHH26fPN3ffto83t5tn5A4fK2ZoiB4MqwffLh58GH3/lFzsxnjX5qhv39DtHfbJwBeu4xU4Hc7MsoXT8mI35M3I/7+1WMy4vd32yekDX/cDw/4PshH+MOvNw/+3Pzd3iaf8Q+/NvDt3cODgA+ZI9ha/Y10YS4Pf/C0/PTe6ibeAB7wpjn0y0gK3iaih+9n2gEe8IY5DMpIBd768PnY4QdTawEPeLMcJmUA3h6ekp/YW8O59IAHvFGOYQAe8P5yzAb/7Qf97/Vxw1NHTQFeBf/+L4+nwt8snu6+fKeaxaHZVMA7LGMMv6ZiD3/7p//8qD3kafjbPxJz5RJDs6nW8PRhkoAfwnNG/Jen5GcS/P0r8h9vpLt19OGH8oC3LEMNf/tkodgjJ4Zv/7Ny079uUwHvrgw1/A2ZavVed8j7XLmzhR8dEA94Kfz9T+RjWrm0Bry/HCl+j78hE7aUy3rdptrBj8+AAXh/8PJpWoDXy5EivN73QO2m9vKAtyzDL/zuvdamH+2m2sAz5zoCvMsYL+r7z/jz5bOz/YPnS3Kruwv4WcoYubke8X1sXlw1P82N7aez4V3AZw5/ebK7e3PR3Lh7+/PyZX932YR29rVxPV3U0/9rCeEanszDbzf7XZ7utp+J9Ob5RXOnu0tC/z3eDXnjYcKezhAj3h/8/Suycrc/CKcb4iQ2J9Rd/aZK4KtKloJzGkvA+4Nvv87t/+4+1DcnZPhP+4wXwVcVYafkAS/P4Rd+MOLb1fi7d1dkrf5kN22tngfforc3xSkAP35ACf/t+/1hspPgB5/xstBv6hi+ohfwgzt0Ct4JiwGvgP/D18OemgnweqHf1CF8NVIn0T8AeHmOMXxFRQf/14n7413Dd/LXLPo+ukepFNwzlAN+CC9Y1GsfRulzzh2JHl7whOPjgJfn0IBvRvzUiRhu59yROMKL3LvfAF6ewy+84zl3Kx34w6+GKfjXogC8Ap4s6qVnuhLDO55zt9KCbzflAF6eQw1vFoFW7ip5UysqheDiM4BPCf4gr4Bv5AEvz5Ep/KoCvDzHFKf54St1U/uVANFlxgA/iERGvAb8rpMHPC9HvvDHMS+8riDgk4Ov9HqrlQc8N0dy8EReE76VBzw3R9bw5IniK8gCPjX4SqepgJfn0IAn22ynTsSYGb55KuD5OdTwZA9bu7MlQXjZtcIBP4SvqSDwX/RPgBMGvtJq6hFevDcH8EN4dsTrnwYlDPzKBL6W7MYDvBz+hox4+YVoYoYXywNeDt9+xusu7/3DV2utpnbwQnnAy+FjW6s3gd+v2gGel0MD3igihBfJAz4p+KrdXK/VW8fvcnx5wGcPz5cHPOAlOeSRDbx1+IYnhLrwg612PHnAFwDPkwd8CfAcecCnBk/kNXqL3j8D+NEDacG3fFPgWXnAlwHPyAM+T3hmTzzgqQeSgj/YTYMfywO+FPiRPOCTg2/klU3lzrmi5AGfDHznBnjrMnKFr/kphvKALwh+KA/4VOB7NMBbl5EpfK1cPwR8ivCr9WT4PgvgE4EfrpgB3raMPOFr5Q6+FeAjgB9P+LpmH2miGtxe857QRy1IwSaShiSHdjjI4aWMKOC13uPUxhebEX/MhBGfHXwtSEGnAnyC8MOrjbKhgtc8/A7wMcCP9qkC3rKMHOFrQYpRNsAXCK91/iTARwA/cge8bRkZwteCFOOEgE8R/losrwevPgE24AEvDMDHA1+LUqhSsgH4LOGFVzCTlMHkUD4D8FbwjNG1cPUO8Fo5soM/zq7V6S2VPOAzhVfJAz5GeIE84PVy5AbfHUeh11tyecBnCy+XB3wy8P2BU4CX50gYnidvDC+VB/zc8KyOM3iZPOBTgR8cIgt4eY6U4Vn5KfASecBnDS+WB3wi8MOTIZj0lkge8JHCj+UBr58D8Pz85jkEAXj/8NRpb8x6iy8P+JnhOSyHplLygDfIAXjhSwA+EXj6BGemvcWTB3y08EN5wJvkALzkRQCfAvzoVJbmvcW+CuCLgGdfBvDxwvfygDfKkQI8b72LgR+ftHhKb41fCPCFwI9fCfARwx/lAW+WIxd45vT003qLfinAFwNPvxbgY4Zv5V3BUy8G+Ojh2QuRAF6eA/DjGLwa4GeF57lTTV07hR+8HuBjh+dccsqit7oXBDzgTXO4qMNrGYDnxPEVAe8C/nz57Iz8u/24XL64au4e7tvDcy8uaNVblYMcLurwWEYw+M2Lq+aH3HjZvAlOt5/O+t/Jy5wB/vCagHcAf3myu3tzcbxzevf25+VLcnPZhDxlpXzRWrs87VC/KIIKMfzpbvv5AN8M+s3zC/JQG/L3p3LE17wDpm2HSeUgh4s6fJURDr4f8ecvD/4nOvBcd8DblhEMvvuM3348PaDrjfiZ4MnrAt4BfLtWf/fu6px8qp+QtfrjgLeEJ6t2rLx9b1WAdwIvCWmZs8GvKsCXCa9xfmvAzwfffoln5AEvzwF4SeyU8oCfDf6w1c4PvHLMA94bvKDrx/CMvKPeUsgDPld4hTzg54Lv9s/4gpfLAx7wshx2dQCe09R+j+xI3l1vyeQBPxP8YE+8N3iZPOBzhpfIA35++JG8094SygPeF7yoy9umUnOuPMIrytDLMbkOwDNNpSfbUfKOe0u1OUEnx9Q6AM80NRy8xnwQdY6JdQB+3NTx7NqhvPPe4lYC+PzhVcduauWYVAfgx01l5tMP5D30FqcWwM8Bzx5H4RdeePZ0kxwT6gD8qKmcA2h6eS+9xVQDeE/w0i0n4eGZegA/AzzviLle3g/8uCDAlwI/qgjwscB38r7g6ZIAHx6e7+4fnqoJ8NHAH+X9wQ+LAnxJ8IOqAB8cXuR+lPcJ35cF+LLgu7oA7wde6C6Db+X9wh8rA3xoeIl7EPhVVdnn0KwD8IOQwe/lfcO3xQG+QHhSHeADw9fypq6DwDflAb5I+FUF+LDwtaqp6zDwh1U8uxyAZ8uMHv5aeeoEwE+AFy/plU1dh4JXyQPeIXytbmo4eIU84MPC8052OQpnK2ZSecDnCy9dxQO8O/hap6k7pbzLr2JiecBnDS+WB3xoeOXC3u3GF9HiHvDO4GutpoaGF1ULeElcj6N9pGIe30fNf5hNsdZ6ojyHSVS8gg1zOChDL0cU8IL3p2TAaw0TxZD3sJ2dUzFGfHh4hbyPHSzsJz3gZ4CXy/vZszYuGvCO4GvNprYpZPKedqmOBj3gZ4GXyXvbl07VDfh54CXy/iZRDAc94N3AH+fazdNb2jn60gFvDC8b8Aa9JRzyXqdNdYMe8HPBC+U9z5erHORwliJ5+G5WtUlvCeR9T5RsBz3g54MXyPufIUvoAT8jPF8+xNToysXVKgE/OH7GsLd48kHmxGtMwAb8qEyX8Dz5MPA7JT3gR2Wy/TU4YM64t1j5UPCqUQ94ukzpgJ/QW4x8OHg5PeDpMl3DMxESXkYPeLpM6ZJ+Sm+Nh3xYeDE94OkyncOP5UPDi+gBT5fpHn4kHx6eTw94ukymi6iTYEzrLUp+DngePeDpMn3AU/LzwBN6ummAp8v0Ar9a9/Rzwa9Gwx7wVJlyd4ve6uRnhKfoAU+V6Q2+G/Szwg/oAU+VOYYfnd/MqrfWeil8bwSqNE+SCHiTpkp7az/oZ4dftfaAp8r0Cr8f9DHAr9iV/EBlpAI/PpOl9TBZryOBb3Io7AFv0lR1b6lPlhJwR4/MvmR45tS1Lj4Y1yr6sHv4hPZFwSvcHa0RKeTDwq/29hx8wJs0Va+35IM+ODwJFr9ceM45yp19B5LRzwJPoqL0i4XnnZve4ZdfMf1s8Pvo8EuF516TwOlWj7XAfl54EvuhXyg8/1oUrjd3cennh98Hf5XPsozo4QXXIHG/nZNDHwk8SVFJ9QEvaak6BbPEjwh+H0L8fOBV7r72bND0scGTqCrO4M8PXnixKW+7tNaDcR8j/CEq6g2QHbz4ImNe92Ue8SOGP0ZVcRcB6hxRw0suLud9J/Z6LfqSZ1SGb/guheIdAHiDFKudEj8e+C7474Ck4GVXkww3bUU69COE74J6B6QEL72KaOD5SmsBf8zwXXC3/sULL3WfZ6Iay58EPDdHnPBK95l6ax/rQQDeNbzCfU74YayZ8FBH4vDny2dngxvd3aTh2Rzqt4JxHWnDb15cNT/dje4uH17lHi88m4N9KzgKq6YEg7882d29uehudHeXTbguopAweZN4L0YMf7rbfr7obnR3Sajfn0wkNOKtc6S9qBeNeMBnDm/0GR9rb82TI234djX+7t2Vzlp9rL01T47E4SXho6mRoEXbFMDHnwPwJk2NBC3apgA+/hyAN2lqJGjRNgXw8ecAvElTI0GLtimAjz8H4E2aGglatE0BfPw5AG/S1EjQom0K4OPPAXiTpkaCFm1TAB9/DsCbNDUStGibEgU8Ew5m4bmYyBdJjkjKkAfg3eeIpAx5AN59jkjKkAfg3eeIpAx5uIFHJBeALzQAX2gAvtAAfKHhAn54rMXU6I/Pmhrbj8tld6zP1DhfLk8tU+wPOrMvw75L5eEAnjq6amqO5XNb+M3LpsNOrXPYvwM3S1v47SfP6jsn8NTxlNNi+0t/KK5VKafWKazhf/v7f23h797+vHxpmUMVLuCHR1BPDSfwG/veOrcdrs0bZ2MLv2kWfw7ew9KIY8S7gT93MkosP6EvyZkjrD/kG3wHOWQRyWe8i2XGx1PLDG1v26+a2Y/4EyefWtKIZa3eHv7cxVBrktgvNuwH67mThYY08D2+0AB8oQH4QgPwhQbgC40C4L/98KH5mbuK2KIM+LlLiDDyh799snjwzx8+3P/0t8Xi8U3zp32o9DdD/vCHRf39q0e7b98/2o//9493Xx7NXdbMUQ78T6935M/tH183P7vbHwsf8kXCP1ksFt+9nruueaNI+NJHO4kS4cln/M3Dr3PXNW8UAH//ql2r7+GbZX3pS/oS4BG8AHyhAfhCA/CFBuALDcAXGoAvNABfaAC+0AB8ofF/UmH7uS3aECwAAAAASUVORK5CYII=" alt="plot of chunk unnamed-chunk-1"/> </p>
+
+<pre><code class="r">
+# function that calculates residual sum of squares
+ssq=function(parms){
+
+  # inital concentration
+  cinit=c(A=1,B=0,C=0)
+  # time points for which conc is reported
+  # include the points where data is available
+  t=c(seq(0,5,0.1),df$time)
+  t=sort(unique(t))
+  # parameters from the parameter estimation routine
+  k1=parms[1]
+  k2=parms[2]
+  # solve ODE for a given set of parameters
+  out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2))
+
+  # Filter data that contains time points where data is available
+  outdf=data.frame(out)
+  outdf=outdf[outdf$time %in% df$time,]
+  # Evaluate predicted vs experimental residual
+  preddf=melt(outdf,id.var=&quot;time&quot;,variable.name=&quot;species&quot;,value.name=&quot;conc&quot;)
+  expdf=melt(df,id.var=&quot;time&quot;,variable.name=&quot;species&quot;,value.name=&quot;conc&quot;)
+  ssqres=preddf$conc-expdf$conc
+
+  # return predicted vs experimental residual
+  return(ssqres)
+
+}
+
+# parameter fitting using levenberg marquart algorithm
+# initial guess for parameters
+parms=c(k1=0.5,k2=0.5)
+# fitting
+fitval=nls.lm(par=parms,fn=ssq)
+
+# Summary of fit
+summary(fitval)
+</code></pre>
+
+<pre><code>
+Parameters:
+   Estimate Std. Error t value Pr(&gt;|t|)    
+k1   2.0191     0.0487    41.5   &lt;2e-16 ***
+k2   0.9930     0.0178    55.8   &lt;2e-16 ***
+---
+Signif. codes:  0 &#39;***&#39; 0.001 &#39;**&#39; 0.01 &#39;*&#39; 0.05 &#39;.&#39; 0.1 &#39; &#39; 1
+
+Residual standard error: 0.0212 on 58 degrees of freedom
+Number of iterations to termination: 7 
+Reason for termination: Relative error in the sum of squares is at most `ftol&#39;. 
+</code></pre>
+
+<pre><code class="r"># Estimated parameter
+parest=as.list(coef(fitval))
+parest
+</code></pre>
+
+<pre><code>$k1
+[1] 2.019
+
+$k2
+[1] 0.993
+</code></pre>
+
+<pre><code class="r"># degrees of freedom: # data points - # parameters
+dof=3*nrow(df)-2
+dof
+</code></pre>
+
+<pre><code>[1] 58
+</code></pre>
+
+<pre><code class="r"># mean error
+ms=sqrt(deviance(fitval)/dof)
+ms
+</code></pre>
+
+<pre><code>[1] 0.0212
+</code></pre>
+
+<pre><code class="r"># variance Covariance Matrix
+S=vcov(fitval)
+S
+</code></pre>
+
+<pre><code>           k1         k2
+k1  0.0023685 -0.0003606
+k2 -0.0003606  0.0003165
+</code></pre>
+
+<pre><code class="r">
+# plot of predicted vs experimental data
+
+# simulated predicted profile at estimated parameter values
+cinit=c(A=1,B=0,C=0)
+t=seq(0,5,0.2)
+parms=as.list(parest)
+out=ode(y=cinit,times=t,func=rxnrate,parms=parms)
+outdf=data.frame(out)
+names(outdf)=c(&quot;time&quot;,&quot;ca_pred&quot;,&quot;cb_pred&quot;,&quot;cc_pred&quot;)
+
+# Overlay predicted profile with experimental data
+tmppred=melt(outdf,id.var=c(&quot;time&quot;),variable.name=&quot;species&quot;,value.name=&quot;conc&quot;)
+tmpexp=melt(df,id.var=c(&quot;time&quot;),variable.name=&quot;species&quot;,value.name=&quot;conc&quot;)
+p=ggplot(data=tmppred,aes(x=time,y=conc,color=species,linetype=species))+geom_line()
+p=p+geom_line(data=tmpexp,aes(x=time,y=conc,color=species,linetype=species))
+p=p+geom_point(data=tmpexp,aes(x=time,y=conc,color=species))
+p=p+scale_linetype_manual(values=c(0,1,0,1,0,1))
+p=p+scale_color_manual(values=rep(c(&quot;red&quot;,&quot;blue&quot;,&quot;green&quot;),each=2))+theme_bw()
+print(p)
+</code></pre>
+
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAAolBMVEUAAAAAADoAAGYAAP8AOjoAOpAAZrYA/wA6AAA6ADo6AGY6OgA6Ojo6OpA6ZmY6ZrY6kNtmAABmADpmAGZmOjpmOpBmZgBmZjpmkNtmtv9/f3+QOgCQOjqQOmaQZpCQkGaQkLaQtpCQ2/+2ZgC2Zjq225C22/+2///MzMzbkDrbtmbb25Db///l5eX6+vr/AAD/tmb/25D/27b//7b//9v///+qozkeAAAACXBIWXMAAAsSAAALEgHS3X78AAAZMklEQVR4nO3dj3+jOHoGcE/cGea2e7eZm51eM5lpr/FeU19j49Tm///XigAD+i0hiRfB83xms7GTF8t8I4yNELsK2WR21A1AaAL4jQbwGw3gNxrAbzSA32gAv9EAfqMB/EYD+I0G8BsN4DcawG80k+C/I4kTm1nONHjh9sny+2fLz6nryRsg1gN+nnryBgCepp68AYCnqSdvAOBp6skbAHiaevIGAJ6mnrwBgKepJ28A4GnqyRsAeJp68gYAnqaevAGAp6knbwDgaerJGwB4mnryBgCepp68AYCnqSdvAOBp6skbsCz469fX5v+H3YeX+9c2gE9cTwp/2T008JePb/0/TbuWtt5868kbsCT4299vPxv44yPr++1XVlDn+4lPeUKihnZTf4d/Yt+1X7ufiO0qLQ+CHu9Zvwx4vser2rVl+KKI0IBFwju8xm8YviiY/Im7x78By4O/fntz2KsHvHiHXwOkilzex28YXtrUT4CXS3KBP1nk1wwv1ftt6osu/L2An6d+3gYU43T3CL8P+Hnq52iAYG2uB/w89WkboOXW1wN+nvpEDRh6+Fm431YP+HnqozdA7OTSa7ylHvDz1MdsgHK7Dnh11gKvfS3Hpl6dNcAbd+Dy+shWH6FdZbltePte+zrhy9Imv2J4+SDNlAYAnqZ+4gKGjr5N+FrecphmhfD81n2j8HW7twUvvaQDXp01wSt35ACvTj7wlvfhut13wKuTDbwO9tz9UPcLloM0gNckB/j2BxPgWQng1VkMvG5TP3D7b+oBr89y4FWZACfUr3hTb5bPGV5xPD1+AwBPU69dwH0TD/g2G4EfXs8B32YL8PoRNAkaAHiaetXQqVkbAHiaemEBlqFTCRoAeJp6bgExTnoMrQf8PPWjBRg+sk3YAMDT1I8+sqVpAOBp6rsF6AfOAb6NavIjk3we8IbxkoBvs0Z44zhZwLdZH7yWvf0B4NusDl7b25UnPcZvAOBJ6otCuwDAc1kNPFM1nxCBTf04a4Hvz32ifq0B/Lz1/Qg6wLtlLfD9ezjAu2Ul8MN7OMC7RTltuUF+kfApz30LrQd8unrurTvg3bICeP4jG8C7JXv41Ge7htYvFv7M58S+lOfJOU0vnVRfUDfAVr9YeOF2Zj0+/WnOofWAT1CvOhAHeLfkDD/L+e2h9YCPXq8+AAt4t2QLrxtvAXi3qC84qJcngu+Z7/Xa8RaAd0se8MPsFV29flQd4N2SJ7xhNCXg3ZIHvLCpN42iBbxbMoHn6pm7diwt4N2SIXzrjr36PhuBN81VNksD/OoBH6f+7o1N/ZAtwFtnKwO8Y/KCt08zB3jHZAVfOJxCBXi3qOH18pTwLmdAAd4xGcE7nQEFeMfkA+82MSXgHZMNPH+QhqABE+sBH1YvHpadvQFT6wEfVO88FS3gHQP4xPWAD6kfPp8FvJz1wntMPg14x+QAPz4eA3g5MeG18gTw3HE4wMtZKbzfdPOAd8zi4YXj7oCXs0p4cbwF4OUAPkUDQusBP6ne+8oigHfMsuHlgXWAl7M++AmXlAG8Y5YMrxpIC3g5a4OfdC0hwDtGB6+TB7xnPeB969VnTABezrrgJ149DPCOWSq87gwpwMsBfMwGxKoHvFf95OsFAt4xy4TXT3kBeDnrgQ+4UCTgHQP4xPWAd68PuUIo4B2jhdfIz7DeTO6AV0QPf9h9eGH/P+7qPNY3u9uKdtHDG+e0ArwiWvjLx7f6X/v97T/ebj9ehp8tEl4/tRHgFdHCHx+r69fX9vs/Xqrr7593n5qCOt9PupTanyRNwf6rQ/Po8UMK/1Tdfrbwl1r88vDK7mqztB7f9nRs6n3i0uP/6Lbyl8fuZwuD118M2vHxAT/K8Bp//fbWoi+1xwN+Qix79Qz9/bf25u7e4RcGXwSvd8A7BvCJ6wHvUF+Er3fAO0YPr5ZPud6KwPrgBqSoB7y9HvDTkjt8EVgf3IAk9YC31gN+YjKHLwLrgxuQph7wlvoisD64AYnqAW+pB/zkZA1/PyoDeP8APsYCAD8rfH8YFvD+iQ2vlE+z3obD74D3D+BjLADwM8KPxtsA3j+Aj7EAwM8HPx5gB3j/5ArPDawEvH8AH2MBgJ8Lnh9JDXj/AD7GAgA/E7xw6gTg/RMdXiUPeM96wKvrxXOlAO8fwMdYAOCrOeClkyMB7x/Ax1gA4KsZ4OWzoQHvn/zgFWfBA94/+cIXOB4fkuzgB/cCQ68CEh9eIQ94z3rAD/Ud87B9x6Y+KLnA3zu4coIjwPsnM3j1xFaA908u8JWhwwN+QrKBb6KZyQ7w/gF8jAUAvokkH2u96aauBLx/AB9jAYBvkgpeO1ct4P0D+BgLAHyTRPD6i1AA3j+Aj7GAzcCf+ZyE2+XZL2K9MkVgfXADZqx3g3//04v9l3RJ0eOlLh+lwxguN4Qe759s4E2XmVo5/O15t/tUXT786+7jWzd9/JHdw3p8c7P9Bb8APsYC0sLXwP/33zX8y/XL4/Hh9frlqfn+qb6/vdn+gl9ygTdeV27l8KxbP7yyqwIdPh2ai8Edm2uH1N6H/tpwD6+WxxAC+BgLSP4az3r57on1+Ia87vG350fW4+9XiKt/wfIYQjKBN7qvHZ69hH98uzz8W/21/Z69xn98q+Hbm92dXkkCL8oD3rNetVd/8d2Ym5MHvNkd8BMC+BgLwPv4LpHhLe5rh/8fOZblOQTwMRaQGF76/Y3AF6nXO+AdA/jE9dnCC/Jh661Ivt4B7xjAJ64HfNXs2QFeCOAd6h0ePxP4shTvCcji4YvAepfHzwO+LHt5wNvrXR4/R/jmMM3BfwDGPUuHLwLrnR4/D3huU//+57fq+JffqtuPiePuEsHz8oD3rLfv3F0e2f/rbv9hnfBFYL3b42cIz3r8hQ23W2mPB3yleTtXd/aHf3zeffib58ibewAfYwF4H9+m3guJA18E1js+PuAdYzuFavTGo8mE9cbNfAJ4IeuFb+c6uh+H3zp8NgMxGPtYHvCe9Zn2+CZh8BU3uRnghawYnqUfcQV4IYA31js/PuAdMw/8MMQS8EIAb6p3f3zAO8YFfiwPeM/6OAMx/jB9jL9c+NFgesCzeA/EAHz28CUX5UCM9798qb8//vrw2s6Zcdg9/C1L+PHZM5uHb2MciPH+y2t1eDp+qi71n8Lxqf56/SvgLT/PBH58jzQQ4/23qr7z+MjuY3NmPGW6qedOlwO8EMVAjK7HP3Z/EnWPN4/RAHyMBSxgIEb3Gt+8uu92T/XXD7/SwI/kAe9ZP+F9PNvU+2Sh8PyJ0YAXMr7n/TN7Uf/wX4uE3++NywO8cHv+T+4u9avD0Tpi1xd+vzfLS/XCTAhbh08/EKN968d2EI0BfOL62Xv87ZntI15sc6Yl39QDfpw5NvXXLzuH2TFT79yJc94AXkjSo3OH+4dCh+bjocPoZB03+EFeft78JgDwwm1K+Hp7327y2w+A+ptyuybACy/6Qr00yRXghaTYq9912/rjY3X9yrb5198/tx8JtTfZz7+fnFIKt/dN7t/p6wq3xa84s8OP5kA+PlW3n0yaTaV5fOpvskzr8V0fb/lNm3r0eDX8sNK6e3w/rRtHgP/a79f1XZylOewz3AyC76Lf1MvTGQKeZdRdosNXh8f7d/2LOjvac3wKf40X39GNbgNeuO0AzwZi/G9zYGZaxE19/xrf7sZfv92vghK8Vy+9lR/u4OoV85cCvsl4U98OxGgOxVqWrUnCz+oHefa8FR/h9HcBXrjtOBDj/d+r7vC7f+aCV350d78T8MJtx4EY/1J/idTj2UUu7Fe5iAV/33iN61UzVQNeyH0gxq+fI73G357ZhuMY57P6Mbz2s/q9WA/4iuADnO4jm69Rjs45wTc/Abxwm+DoXJoebzo4t+fqldckALyQ+B/Zxn2Nv8ubD8ruAU8P75a48NV+qFdfhGTr8NlMhdKnhd/bnvfwhwH4JnmPuauc4U+9POCbzL9XH3PMXdXB7x2edyevudwU4IUse8xd5QPPzWsnBvBClj3mrmrh9y7Pu53eDPBdst+rr+XZRtwVXndhQcA38R2IQXS2bBNX+Fp+D/g+cQZikMI3TXV63nv9lUQ3D7/nohyIsZwZMZp4wNc7eLqfbx6+jXEgxnJmxGjb2nx1fN6Av2fCQIzlzIjRxAe+0H60C3ghioEYy5kRo8m+eScPeM/6CQMxFjQjRsU6vBe8Th7wQpY9I0blA79vd+3U8oAXsuwZMSoP+P49qlIe8EIW/8ndvv3UFvCe9bnDM0U3+OE9qkp+6/DZDcTwgB/exCvktw4v/f6K4Ecf3sjygBcCePPju9YD3jGO8K1g6QsvywNeyHrg+Y/pRXnACwG88fGd6wHvmOjw4nE5QR7wTYpCvCcg6eFreW94QR7wLN24xNE9IUkIf8cDfAr45po0R4fz3TRZBLxiBAY3LRbgm4w39e0IHPZleTNihMDzE6IBXkg3AmfqNCgsgI+xAMIROJeJU2Kkgx/gSsvzLhT1gB+iHYHzenQ4+0WTpcJz+3eAF7LovfpQ+NECAC9kHfCF5XMAwEu/nwk8d0lxOYBf1UCM8fvwifD9MrYOL/3+KuAL85g9FsALIYM/8zmd5exH35eKn/cp1PXSQgzR1rsmdAGx6xcLL9xO1uP9zr3TBj1eTip47ijL2SBfqOu55QBeCOC1j88F8HKWDu86o4YxK4H3HIhhPEuaHL7Q1HNLAjyL70AMEnhhIAXg/epF+IKLaiCG70woM8Hr9+vt8E7z5JmTPXwb00AM35lQqOELTT2/MMALkQdi+M6EAvgoC6AfiOE7EwoxfKGrNy1NzgbhxYEYvjOhJIKXx8Wr5QGvrvd/H+990cl54MsyCN54aROXesvPc4afOBPKPPClBv6+l2pfb2b5LcNPDeBjLADwbaSTHoPhzfJrh89mBM5e3rlT7d31n0e4rDeT/MrhkyQFPH82RAV4wPPxgjfJA94/hJv64cQpt/Wmlwe8f5L0ePH3m+ctyQNeWw/4Ub1BHvD+oYMfnSLrut508oD3D+BjLADwLJoJCjn4opgCr5MHvH/mgx/Lj4aPyfVyhvWmlge8fzKDV8sD3j808JU/vGk+e8D7Jz68dtpx/kVeW6/IuTJfyQDw/pkRnu/y2npFxvCq5QPeP7nAGy9hAXj/EMHzE5z5rTf5AQDvnwzh5UcAvH/mhB/JB8G3DzE6BAh4/9DAC1NZToEfH/QHvH+iw5suHxYLnj0I4MMyK3wvHwrfygfUCwG8W+jh+YcBvH9I4MXZqgEv3Aa8UK95HMD7Z174Tj4G/PiBAO8fCnjpugSAF24DXqwfB3v1IZkZvpGPBF/hfXxAYsMrB0rw8PIVaKaut/uDAd4/gI+xAMBb4Wv5ePD3RwO8f+aHV1xdcPp62wfWBzcgTT0p/GH3oZk26fbcTJ942HW3Fe0CfOR6SvjLx7dLM10mm0jt8MTNmhYCr7qeaMB62wfWBzcgST0l/PGxun69X8zu+HT9/fOumUqPzbH0/aTNXv+jLoX1N7xif8AMQwr/VN1+dvB1p788vLK72oT0+EIxQ0JIh9kH1gc3IEX9Qnr8obt8aT916pLg2SMC3j/21/jbM+voDN2hx6vPbePgFXNjhK23PeAnxLJXf/32dmimQ67aadDV7fKAL1STogSutz3g/RP5fTzgo9SvE16WD11vtkmOAS9nZvjmTXx0+JNFHvByiOD5KU7D15tZHvByKOCbSY3H8hHWm1Ee8HLmhe8+rk0DL06r6F4PeMcEw0ff1Ldn1ujkAS+HBF7cvYuy3gDvlbjwmjU/Czw29V6ZFb4/IpsEXr+DB3g5a4LXygNeDg28IB9tvdleaqwLCG1ApPqVwY/G3iSCdzg6aFlAaAPi1AOeq9cG8P4hgq+iv483tQDwcmaE50ZZpoK3ncJlX0BoA2LUA35crw9Xb5h8yW0BoQ2IUJ8bvPntFD+uOvJn9aPoJlF2XkBoA8LrVwUvjKdPB6+5QobHAkIbEFy/ZvixfOz1Jl3i1FIPeMcsHV5x/fr27lQXtQG8Oux5y2dODfLx17t8pcvKdPwO8G5ZPjzfuQEvhxB+kE+y3keNwaZezlzwCvfE8KPWYOdODiV8L59ovffNAbycmPCmgRAU8H2DAC9nJniley+fbL3v3eoB75hs4LsmAV4OLXwnn3C9713qAe8Yb/iCDL5pFODlzAIvXEp2nHJUr03Qeq/fvANezvrh63YBXs5c8NpllVX69W47fx7wjvGGNywrIXz/Ce3ZJg94t6jh9e7G510m+8i2PyZztvV5wLslP3jDiXX6BYQ2YHr9JuBr+eSb+uaW/wJCGzC5fh3wheV5J4MX6g3ygHdLXPiqnGm96+UB75ZM4fXygHeLF3xhf96K6W25RFvvGHrVZ1vwukYC3i2x4U/ibEhCYq53ZSsB7xYlvN7d4XkL858JibreVZt7wLvl+5nPiX3Zn5Up1HcL9Qze4Re19V6RW+q5gOT1i4UXbof3+Pk29Sx9U++9Hz3eLR7wzfEZl+c9J/wdvP9QF/BuSQFvkk+w3pvWAt4zK4C/y09fQGgDTPX5w7dH4t2et14+yXof790D3i1p4PXyidb7QA94tySC18onW+/Y1PtFBW9wXy78nR7wbkkFr5NPut73Ea5YCHg+90GW7s9bLZ94ve+Dr18GeD7+8Gr51Os9+PplgOfSj6pePPzJPBYT8F0Swivl08NbhuECvk1KeJX8HPDG7T3g27jBD+fP+D1vWX4e+ICL2gB+nKnwjTx3mHYm+OaojfCTvdMJ9huGN7l7P+9SGJEzJzz/RLo7AN8mNbx4KcpZ4Tl6wHNxgR+fIeu/3kqaTX27ZR/TY1M/Tnp4otf4e8QtPuDbOMBzp8RPWm8Um/pROHvAt5kFfiRPAV+NN/mAbzMP/CBPBD90e8C3scPzk59MXW93eTL4quv2Z/5m9AYAXkgZWO/4+OYF1N3+PPpeJQ/4LsJsR9PhysB6t8e3LqDnBjyXk7nDh8A1b+vI4c+jczBSNGAt8OL0ZkFwNf0S4FUf5EdrAOCVKZcBX+ntAd9GnMgyFM54RqVDfTT4SvMqD/gm0tS14T3WLD8nPItkD/gmCeDNnX5u+Ers+IBnKaJv6lm9acoMAniWBr/xt3wOINTLrxVrgRcTpccbZkshgmdpe76pfj/+AGgoMT5+lvCKqaojwWs3+MTw0vAN6TdWC2/u8JE29c23Q7cf/REQwg+bei2/fVMv3ZEjvGpu+rjvw1v78Zb/1N3vVh/cAFO9mt9cL1dkCK+8JkH0D2DKUoY37APMCN9G5F8/vPpaFAk+eeOYFwffZnjxt9Tnv6nXXIMkzUeug/SCNvWKmPf9lPW5weuuPZPus/YWXzeCx3Ugh0cD1H9dTkf3DP6Zw2uvOZT2IIvuHf5wdzx4+aFK8WpKKtzxCB4Vf97w+mtNzTBkTYE/C7y0yVF2a7EB4rv/rOEN1xibaayiqGLf1PM7CS4NkDv8FHj+1+UZOXKBt7jPOUi1vMelvv89w9iv0tIAz029OhnDm9xnH53MUo6jqbPDd78R+xMoWz0p/GH34WX0TX9T0a692Z0EnqvX/RVYN/Up4BV/iUuCv3x8q//13/Q3Ve3am93p4bv/tyu81EVVat3UezWguj++pZ4S/vhYXb++9t/0N3d1vp/4FKcswta47ec0EVtCCv9U3X6+9t/0N1k0FxXWZiE9Xterl/cEltjjVe1a2nrT1ScbyTGuVz1ITvA+r/G5wM/RAOVmJSf4djf++u3NvlcvDrGTAnjPx8/jfbw0qFbKluCz39SbAvjE9XnAr2lTr9n9A3yb1e7c6d7wJW6A9KCAn6eeGF5+VMDPU592U289zztbeMuAspzgEzSAsa5zU28dSgj4de7cAT5wU58rPDb1eDunCeA96wE/Tz15AwBPU0/eAMDT1JM3APA09eQNADxNPXkDAE9TT94AwNPUkzcA8DT15A0APE09eQMAT1NP3gDA09STNwDwNPXkDQA8TT15AwBPU0/eAMDT1JM3APA09eQNADxNPXkDsoEXshPv8Ax1PXkDpPrYzHImwcdeCHU9eQOiKBA8JPXzBjzRQ1I/b8Bn8ZDIEgL4jQbwGw3gNxrAbzQR4MfzYU3KMIHelNyed7t+IrZJOex2TyH1zVyAYY8fuAYnJByemwFv0gJ2DyHwl0/1qguBqxcQ9qdXP4MQ+NuP2dWrGPDcnJcTcvv7MFfq5DYE9tgw+Pc//zME/vr7592ngPppiQA/nuV2UoLhL4Hr7RDUY+u/mktI/aXe4IX+5fqHvseHwx/C+0vIi/SRTegd9CJf4wfW+2cBr/GB8LfnwN7CVnrY3lmYW/P4Gfb48L36MPhDcIerlxC4zQjrsIfgDcaE4H38RgP4jQbwGw3gNxrAbzSrhX//5bX+R92K5WbN8NRNWHTWCn/9snv4xy+vtx//Wb9JvjRvlNld+GO4Z63w3ab+9vypev/8qen/h8fqOP/BkKVm7fA/XprDnte/vtT/Qg8qrCkbgv+yoxjwsNRsCB69fZztwLPX+MDDiGvKauFvz+1e/QBfb+uxpe+zWnjEHMBvNIDfaAC/0QB+owH8RgP4jQbwGw3gNxrAbzT/DxiY9hkZGetKAAAAAElFTkSuQmCC" alt="plot of chunk unnamed-chunk-1"/> </p>
+
+<pre><code class="r">
+# Get the 95% confidence region
+
+# Inverse of covariance matrix
+Sinv=solve(S)
+
+# draw the confidence region
+# get points for a circle with radius r
+r=sqrt(qf(0.95,2,58)*2)
+theta=seq(0,2*pi,length.out=100)
+z=cbind(r*cos(theta),r*sin(theta))
+# transform points of circle into points of ellipse using 
+# svd of inverse covariance matrix
+Sinv_svd=svd(Sinv)      # inverse of covariance matrix
+xt=t(Sinv_svd$v)%*%diag(1/sqrt(Sinv_svd$d))%*%t(z) # transform from circle to ellispse
+x=t(xt)
+# translate the ellipse so that center is the estimated parameter value
+x=x+matrix(rep(as.numeric(parest),100),nrow=100,byrow=T)
+
+plot(x[,1],x[,2],type=&quot;l&quot;,xlab=&quot;k1&quot;,ylab=&quot;k2&quot;,lwd=2)
+points(parest$k1,parest$k2,pch=20,col=&quot;blue&quot;,cex=2)
+</code></pre>
+
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAAclBMVEX9/v0AAAAAADkAAGUAAP8AOTkAOY8AZrU5AAA5ADk5AGU5OQA5OY85ZmU5j485j7U5j9plAABlADllAGVlOY9ltf2POQCPOTmPOWWPj9qPtY+P2/21ZgC1/v3ajzna24/a/v39tWX924/9/rX9/tr9/v3lmJ4qAAAAJnRSTlP/////////////////////////////////////////////////AKd6gbwAAAAJcEhZcwAACxIAAAsSAdLdfvwAAA6dSURBVHic7d1tW9vKEYDhY8IBmqYtJIWcFuJibP3/v1hbNuFFtrzSzuzO7DzPh0Oui8NqtXckS46N/+goZH/UngDVCfigAR804IMGfNCADxrwQQM+aMAHDfigAR804IMGfNCADxrwQQM+aMAHDfigAR804IMGfNCADxrwQQM+aMAHDfigAR804IMGfNCADxrwQQM+aMAHDfigAR804IMGfNCADxrwQQM+aMAHDfigAR804IMGfNCADxrwQQM+aMAHDfigAR+0RPiXr79050GFOwu/vlvs+wJ9S50/4td3W3KO+NZKOdWv767+B3xjpT3Gv9xwom8sruqDBnzQzsO/3Jy6qF+Q5TLhNw/3/dfV1fMA/uxfGqpXLvz6x68PX9OHpqpxxActF/71qbt3j/GJjyJUtWz42UNT1YAPmhQ8F3fO4ogPGvBBAz5owAct+5m706/AqQif/px02LKP+M3D7byhxcv9d4lY5Z/q198fZw0t1Dli+I/n+TF+qij873IKnyEIf59HeAm28Pze4GWtAuO7gldCConvB14XJxy+E/giKqHwPcCX5AiDbx++vEMIfOvwtQCat7cNX3ftm7a3DG9g2du1twtvZcUbtTcLb2mxW6Q3Cm9uoZuztwlvcZEbo7cIb3aBW7I3CG95cduhNwdvfmEbsbcG72FRm6C3Be9mQf3bm4L3tJje6S3Be1tI1/aG4B0uomN6M/BeV9ArvRV4p8u3yye9EXiPS/eWR3ob8P7W7VP+6E3Ae1u0Y3mjtwDva8VO5oveALyn5RrPk3x9eEeLdT4/9NXh3axUWm7O97XhnSzThJzQV4Z3sUZTcyFfHT5jA3ZzQF8X3vzyzM3++b4qvPXFycm6fE1440uTm236uvAZozvItHxFeMOrIpVh+nrwdtdEMLvyNeEzxvaTVfpq8DaXQyGj8rXgba6GTibp68FnjOwti/KV4O0thG726IEvkzn5OvDGFqFIxuirwNtaglLZkq8EnzGs4yzR14C3s/elMyQPfNnM0FeAN7LnlbIiD3zxbMhnw7/cLC4ep3yMuIndrpqJgz4XfvNw33/YZDK8hZ2unQX5XPg9+NP1BPi0ibVdfXqJI37b8s+vwE+punz2Y/z6rv9Q4WXip0nj/lpl+dJX9cD/ru5BD3y9qspL3M7tGp7pjw6N+4cqygtd3HWrq+eUoYH/WL2DXuZ2LvkJHOA/VU2+7BGP+7BK9AK3c58f4xevHdka8MPqyJe9qgf+aDXki8LjfqIK8lLwSRd3wJ+qvDxHvI2Ky5eEx32k0pd4wFupsDzwdioqn/3M3d3hpv38P8vifq6S8tlH/O51V2lDA3+2gvL5p/r198e0oYE/Xzn5go/xwCdUTB54Y5WSB95aheSBN1cZeeDtVUQeeIOVkC8Hj3t6BeSBN5m+PPA2U5cH3mja8sBbTVkeeLPpygNvN1V54A2n+aIc4C2nKA+87dTki8HjPi8teeCtpyQPvPl05IG3n4o88A7SkAfeQwrywLsI+KiJLx/wPhI/2QPvJGl5nrL1krA88G6SlQfeT6LywDtKUh54TwnKA+8qOXngfQV81KSWEXhnSZ3sgfeWkDzw7pKR5/3x/hKRB95hEvLAewz4qOWvJfAuyz/ZA++zbHngnQZ81DKXE3ivZZ7sgXdbnnxBeOSFAz5qOQsKvONyTvbAew74qM1f0pLwyIs3/2QPvO+Aj9rcNc2Gf7k58UGTwBepFvzm4b7/urp6Pj808ArNXNRc+PWPXx++jg4NvEIzr++KHvHIa1QH/vXjpJMe44FXadaqFr2qB16lWSf7svDIq1QFfsrtHPBKzVjWshd3wOs042Rf9HauQ16p8vATj3jglZq8rgq3c4vXjm4PeJUmn+wLX9UDr9VU+dLwyGtVCT714g54taatbPEjHnmtpp3sgW8n4KM2ZWnLwyOvVkn4w238sSfrgS/ehKXNPuI3D7cThwZerZLw3fr748ShkVcrfWkrPMYDrxfwUUte2xrwyOsFfNRS17YKPPJ6pT5xC3xrAR+1tMWtA4+8YsBHLWlxgW8vy/DIa5ayuCnwLzeLi8ejr66aPzTwignB7146v/s3OGF45PVKWNwE+D3407UovLtD/vLysvYU0pOBP7xZZvnnV0l4Z/KXl43JpzzGr+/611osj7wjNmNoV/CXl77kZeDnbvzc9x3Je4M/v7hJ8Ov/7P67+Uv0iAdeMxn4/nV1y4sTL7GaObRD+dqTmNK5xU071W8e/nFz6jWVc4f2Be/rqr6Tgt/ezU12Pw/vS95bZxY35T7+9Evnc4bugFctH372lhP+F+T1Gl9c4JvNMjzymo0uLvDtZhkeecWAj9rY4taGR14x4KM2srjV4ZFXDPigmYZHXrHTawt80xmHR14r0/Ac8oqdXFsL8MjrZR0eeaVsw3PI63VqaW3AI68W8EEzDo+8WidW1g488joZh+eQ18o6PPJKOYBHXqUTHwA4/kMF4TnklTIPj7xOHuCRV8g+PIe8TkeX1RQ88iq5gEdePgfwyKt0bFGNwXOy18gDPPIKqcC/3Cyul0d/W8ZceOTFO7KmufCbn4/d8nrr/+154tCnAl4+BfjdL7hd3h791dZzAZEXz8URz8lePgX43WP8reRjfIe8QsMVNXdV3/8k8MI5gUdeOg34p/6XH94f2VjqrI78KPKyKcA/bcmX993T9XBj6fMa/izyosnD97dx20t6wdu5/Q8DL9pgPSVu57rVdbe6ErudO/w08pKJw/e/4vjqefXul9kvXpsxv3cTQ14yefiRbWX8bIe8bI7gOdlLpgYvfXHXD4C8WJ6OeE72kn1eS8vwyAsmDv9yc+pjSwTQkBdLGv7w+aPy9/GHMYAXShr+9aJO4+KuHwR5mZwd8ZzspRJ/jD98OJnOY3yHvFSurur3wyAvkT945EVyCI+8RB7hubQXyCU88vl5hUc+M5/wyGfnFB753LzCI5+ZW3jk8/ILj3xWjuGRz8kzPPIZuYbniZz5+YZHfnbu4ZGfl3N45OfmHR75mbmHR35e/uGRn1UD8MjPqQV45GfUBDzyk3P1btmxsZGfVivwyE+sGXjkp9UOPPKTagge+Sm1BI/8hJqCRz69tuCRT64x+C77F2gGycvvq5+wEeRTag+e031SLcIjn1CT8Jzuz9cmPPLn0vj4sZGtZfzs5G0hP1a78MiP1jA88mO1DI/86VQ+Rnxkcxk/O2t7yJ+ocXgu7k/VOjzyxzu6KE3Bc7o/WgR45I8UAh75QccXpDl45D8XBZ5LvE+FgUf+QyfWokV4TvfvCwXPQf9WLHjkXzu1DK3CQ38oHjzyfWrwLzf3m4fFYvh5g9XhucbrTruLfNLk0/3W/5vOJ03mxUGvBr/+8Wvz81Hvs2Uziy5/eu+zT/Xbw31123Wr66lDFyo2vSJ899R/tuzQ3Qh8aPmRPW/4qv53cemDw4eVH9trqVP9/eShixaTXhV+dy+33N7SWb24O7QISD+6wwK3c/1NvNHbuXfFk1eF72/it/dyq+FTd+aWORj9+M5mn+rXd7vna1cXj+82eChtfgUzOSm1lOFnD12lQPRndjQYfBz5c/spBW//4u61GPRndzLaEd8FubM7u4cB4SMc9Of3T+CFGP0R9GVwpjcM3zx9wt5JvBBjl4f7+Pc1fb5P2TWRZ+46Txd3r7VLn7RfUY/4Xa3KJ+2VyDN3/h7jDzVJn7ZLIa/q32rwfJ+4Q8Hh26NP3Zvw8I3RJ+8K8F1L9On7AXxfI/QTdgL4Qy3QT9kD4H/nnn7S9IF/l2/6aXMH/kOO6SdOHPhPOaWfPGvgB3mknz5l4I/kjX7OdIE/miv6WXMF/kRG3xowbOY0gT+ZD/q5cwR+LPP286cH/Hi26TPmBvzZzNpnTQv4hEy+CTRzSsCntbCGnzsZ4JOzZJ8/D+AnZcReYA7AT62+vcjmgZ9RTXupTQM/rzr2glsFfnalL/ZkNwd8Toti+OLbAT63AvgaGwBeIk18pZGBl0oBX/Xv0/i3gZ+UnNRioaneAS9ftpi6+X4r498GflaLmXhlzPebGv828PNbzK3M7Ma/DXxuBs338xr/NvCtBnzQgA8a8EEDPmjABw34oAEfNOCDBnzQgA8a8EEDPmia8GQ5PfijfxkMjxZpcmeHA97IaMDrjRZpcsDrDWd6csDrDWd6csDrDWd6csDrDWd6csXhyUnABw34oAEfNOCDBnzQgA8a8EEDPmjAB00I/uXrr/3Xm8WX7Z/Wd4urZ7HRlotF/3XmYDeLxX3/p8OssiY3GE1qcoe9zlu5wXCnZycDvzqMvr6775ZXz5uH7ZdrqdG6p/tzPzHS+vtj9/K3x+2fDrPKmtxgNKnJHfY6b+UGw43MTgT+6eK/+2P05dvzbuvrH79+H7X5o21+PmbMbbVbxn7/D7PKmtxgNKnJHfY6b+UGw43MTvZUf6A6fBEabXv2ezuBzWo/F6HJfRxNanLdfq+zJ/dxuJHZycL3J+eLx9WVCPxhtN3JK+vA2jzc7r4cZpU7uY+jSU2u2+917uQ+DTcyO/mLu3/9lDriD6P1f8x4KF3f7ddC5oj/NJrU5DqZI/7TcCOzE4bv+ge/vEeqz6P1X+ev7cvN/dtQuY/xg9GkJtcdzs15K/d5uJHZCZ/qt39Z+wvn25yr+k+j7U5/m7+ypV5nlTW5wWhSk+v2e523coPhRmYnCN9vaZF/qzwYbXszejH77Lfs31Vyf7jUyZ3ccDSpyUncxw+HOz07nrkLGvBBAz5owAcN+KABHzTggwZ80IAPGvBBAz5owAcN+KABHzTggwZ80IAPGvBBAz5owAcN+O7tLQc5Lwn3FvBv8KuMt726C/gefvNw9fz7zZohAr6Hf+rfecSpPlabn38/vB8O+FBtHv757/7dK8DHanuqX3KqD9ju4q6/sAc+Vjv13hx4aj7ggwZ80IAPGvBBAz5owAcN+KABHzTggwZ80IAPGvBBAz5owAcN+KABHzTgg/Z/nPbhGH1cTt4AAAAASUVORK5CYII=" alt="plot of chunk unnamed-chunk-1"/> </p>
+
+<pre><code class="r">
+
+# Simulation based estimation of uncertainty
+
+# store original experimental data in a separate dataframe
+dforig=df
+
+# conc profile based on estimated k1 and k2
+cinit=c(A=1,B=0,C=0)
+t=dforig$time
+parms=parest
+out=ode(y=cinit,times=t,func=rxnrate,parms=parms)
+
+outsim=matrix(0,nrow=nrow(dforig),ncol=4)
+outsim[,1]=out[,1]
+
+# number of simulations
+nsim=1000
+
+parsim=matrix(0,nrow=nsim,ncol=2)
+colnames(parsim)=c(&quot;k1&quot;,&quot;k2&quot;)
+
+for (i in 1:nsim){
+
+  # Simulate data set by adding normal random variable with mean 0 and stdev from fit
+
+  outsim[,2:4]=out[,2:4]+matrix(rnorm(3*nrow(dforig)),nrow=nrow(dforig),ncol=3)*ms
+  df=data.frame(outsim)
+  names(df)=c(&quot;time&quot;,&quot;ca&quot;,&quot;cb&quot;,&quot;cc&quot;)
+
+  # get parameter estimate for the simulated dataset
+  parms=as.numeric(parest)
+  fitsim=nls.lm(par=parms,fn=ssq)
+  # store estimated parameters in the ith row
+  parsim[i,]=coef(fitsim)
+
+
+}
+
+# plot the parameter estimates from the 1000 simulations
+plot(parsim[,1],parsim[,2],xlab=&quot;k1&quot;,ylab=&quot;k2&quot;)
+# overlay the 95% ellipse computed previously
+lines(x[,1],x[,2],col=&quot;blue&quot;,lwd=2)
+</code></pre>
+
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAAclBMVEX9/v0AAAAAADkAAGUAAP8AOTkAOY8AZrU5AAA5ADk5AGU5OQA5OY85ZmU5j485j7U5j9plAABlADllAGVlOY9ltf2POQCPOTmPOWWPj9qPtY+P2/21ZgC1/v3ajzna24/a/v39tWX924/9/rX9/tr9/v3lmJ4qAAAAJnRSTlP/////////////////////////////////////////////////AKd6gbwAAAAJcEhZcwAACxIAAAsSAdLdfvwAACAASURBVHic7Z2Llts4coYXjnfsTDbJeDb2bGJv77htvP8rpiURdUPhRoISqar/HM+oSeJCfqgLQEr8S3SZ1F8e3QHXY+TgjcrBG5WDNyoHb1QO3qgcvFE5eKNy8Ebl4I3KwRuVgzcqB29UDt6oHLxROXijcvBG5eCNysEblYM3KgdvVA7eqBy8UTl4o3LwRuXgjcrBG5WDNyoHb1QO3qgcvFE5eKNy8Ebl4I3KwRuVgzcqB29UDt6oHLxROXijcvBG5eCNysEblYM3KgdvVA7eqBy8UTl4o3LwRuXgjcrBG5WDNyoHb1QO3qgcvFE5eKNy8Ebl4I1qC/jgOrJ2BL+hrGtvOXijcvBG5eCNysEblYM3KgdvVA7eqBy8UTl4o3LwPWqucJ5PDr5DIT7T2dzk4Dvk4DdUfWZv6eDXV33ua3fmUVuQgzcqB29UHuONyrN6o3LwRuXgjepc4D1RmKZTgfepwTw5eKNy8K02nzS6nAr8Ayg8rZM5F/j7y8FPr/occvDTq36MhoOFx/jZVT9ET2vAw3LwRuXgjcoY+KcN2cOyBt61yMEblYM3qoOD95C8l44N3pPw3eTgjcrBG9WxwXuM301zwL/++s/hql0P1VbwPz4tv5D5bxl6B39kbbb4H5/ekD/S4o8ZDY7ZK6IJrv7Hp1/+9Tjwx8z/jtkrqikx/vVj7ugd/PF6RXXwrL6vleNd4mP2iur04JdoerSYerT+ZNoM/vVjIam/64g/voUdTVvB//zy+fr/77/8OVj1VM0Af3gjnarN8/i//5P9v7/qqZoA3pjTeBKLn2CuDp6qZwFHxvjO190cTQ6e6gRZ/TSdbqhukoM3qlngH5zcuUblFm9UDt6onh78+5Ie3bEH62nBF4H3jYCnT/Gf7wmcClc6Va/Df/5J/WaL//nlt3VV76GGIas4m0PkObXd1f/4/euqqqerHbyLODP4Dn6/qierK2OrRm4GvxHjz58CPAf4aXl6Z87/BA7hGcBPnp418v2LHPyGqidpE/Wixwb2+n4Hv6HqrioaobTLL5fravCr1O4xfn3VrdJvl7YLzEBvZF0tww1DI+tcOir425Oz1VrGiKwD//bvSdGfFvwojaD555bHXvY/I/pDg1fB3LaNo7iWW9un5zP7o4IvGiO431X9Wd+nZ0N/WPC1Wtch2DwHeyr2JwS/+upvn4M9UZJ/OvDruc9r/hnYnw38ES76U6A/GfijXPHzoz8X+NVXe/4S69nJnwr8eu5R786m8XBu9GcCLy90P7YC+I0TvFP7+xOBz7j3t9EBvrxIWO/SWdGfB3x2iUfsVUdIalhzD2fp1TnRnwZ8fn0nPA2xjIflO91rwJ92cncW8Nq1nfWzRyGut/ilb+dDfxLw5Qs7w+4jGL3Y0T2mzof+NOBrrfS3xFEmTz9URUFnQ38O8LMeeIWD4eY8/rVZ5yJ/CvDVKzpCLbB/Uyyd6kzozwB+4+UkI2Nn8Gfy9+cAv6U0w8uj+g4PSTP0R34I+wTgLxdywyVU7ZrH+KlC8vvUP0nnAL/lEpbL6nsKY2xg6CX0Dn6Trga/pboiMt0X6E0N9WAh7+C36P0g+Lpl4hpt6dj0DR65ub8HF93QP3WMf/0Y3n3d8Xfuur6uzhuVqRx+Tknd9X9qndejtoPfZWY3dRxtBX/5EePLr6HsBX74axOk4fwzyeYLX9MpjIcVV306+rmRYyv4G/BvH84EPsBMTgNf7Pgo+tnkjwV++dnyl7/m76Gasww6WoL8brYyCJaf1a58FTfgsBA79OPL3di88rSpAz21VdTzs+XXX7162efnzlZ/Uyrjlz4vMRwGgdbLkXS/1onJRn+oGL++6h6t/qYUCebigtENQ5ncGvAHXr1/avAso6O70+ehFH7N3aDZ7n6eZkzn9nsL1aqLhgH89klgRH4Bsnutim3CGg5q9JOSu33eSbPlR41Cyt6LPvqW5DV6OSGuHpP8nOncTgs4678/sSzXlJ7KC8nJN7nHOedxPPJHtvgNBl81dTadb9QU0+xvkw5o9BOmc7u9hWor+PqvasQOR4+D4+nIHzmrX3+tCndaYHfEOCBKlWqqn0zHKD8a+ecEL6buAsuS+kXoIq7f6rP3lmvocgkHM/pZ4HdI7mZdJz5zB+q4OdQMuydu9cWCY5E/sMXvAD4hDgJ8iPrtuoorKLSg7w/HW8x5GvDlZysK4PHQAOgl96VIbJ1M4+kPKH8goz8u+OEfriwy4mt1pawuywOwto0zFNKp45C3AJ4cRQO2du8ufV7yfqiN7B8fBaxTRyG/Ffx+b6HC98P0XOlO8GR3ltXT7TjTJ9uCrEJUXuwn23MQ8pstfre3UA0+nV6O8aJLBLzKElbw+WBgrr9Vd1XHIL/d1e/1FqoZX0uordGklE6G8WXJj5UkD3YIX0CP6e3nIQL98WP8BvBLpC7deaXP3S60lz3YYLL+dCsvmwuKxrp0BPKHBb/9G2gJmdoV5rX5dJ3FA3TycEdvKZ7lhiO3Jx5P/gTgaZWluy76w5GVWyyLGctMT7p5Aj5y8ClO8MrWk59wU2tIRwR/vQYa+JI3VbeDh1YLya2Q0IlMD4YOWwAIHHMqKe8M1WAK8lsymVU6IPjbNZgD/vpBew4j88uBwAx8M/mM/2cBIarg6zA5eQc/Dbw+wQO+AicOE60uMdEPOOfDWvMRVoUZAiPv4CV4fs0HYjyrDf+CmRnfny/oCVMP7DOzbbH0X2g638fJe4xnMb5qzD3be8HHKIwYy7Gsjgf7chdr3UyFHpnbHxH8RdVZfMXl50s25FiwVraOwyxebr8Vu/67ERfjo9KbqpJbexj5JwOfZ9bEESe+QeymMZ6CJ04/zQ54elDvjS4eSR5H/ujgS1P0KOndtivgsYQ644oRpmq3apXcHT+q64Aj4Vl2/WHkDw9er7rgbIuLKHXwdAKAZIQjmPOkdTZmH0X+nOCVljD5Kh+hTLqWESGSPTZTD7h1+4wrjwsPIv9s4IuHkEm9nCFqWb6I8Z2t8MY69zl4qlUP3PUcVFpapfN4jpgm/F3nNJrmP4T8scFvX9VAS8dQLcFnj16I9B7u23ca8jj4B5A/KvjtP28HfZDhOksAU5wnsMmkT9yIy9f1SDu4Z6jbjyD/FODLpkiidb4oG8koKLj3mHK6wOrLeibmAqMryw8gfx7wTbrlXWm9rbCantI71VsHYf5d4If7eX/yhwbPbLNco7ZLLMOEUAYvBgfZIW+/sniRLwDKFFC2V3Ngdyd/bPCitn7wEldAirctZIZGwaddIa3m4ew90HtwyihiOYTarWrkujf5ZwDfviPKYzxN8mBEJPDL+MDpvTqLhxhBq29yrk5S7kz+NODHpnYZeCCa6iJuHGxZoAYfkNdIXAWMgpjNEgvXoHAm9yV/WPDNZQ3l8lGbFguzLMTTcQDMwZ8DweTlKXhm5Hxz7vxry8fqGd+R/GnBs8uXbDZGzlsZBwvasrmSFJ/hxC2yCyEq9t7quTYq7kn+LOCLuRT5LMFTg06VQLRnGxm3IGonPbj9yzql2ntRLEzIc74b+ZOAVxKl2As+WSr5HDgq8VSOnPvhVK0EmA2vplgIkSd9L/KnBa/N8cW1R/DJqHEQcJsjI4EGAYI7GXspboux0lQp0Dv4JnjWVBFHlOABo6gyPY0XkWKizxP+QhbR6mJnl+9G/iTgxyZzsiOQrxGfjYad5ulpQEC+z8GjnydBHv7KBtJa3cvZHxd8YeyvGwHCwjFxJ/M1tqCDcZgkg1lspqFhzNWXdSfyZwJfSYoaPaGT9tsGSAPJ5syo0T2wKE7BpxASSQ64UfchfyLwlWlQ1rRM8iSyNIZkaCYhnjsDGBSY8EFDARJArXMkZezUOcDv93v168FDLs6W4zgsZW0ngr9Gthn41ACkBsRRSF8x1uvyme+hreB3/PVq/TZNw3bQBIEDZPW0GnK08iGQIUDzQawDaqTEg/gbmntC8Dv+Xn3jNo06AmBwMByIMCpQtFka5H+QDrAFH2r+tM6UAK4AzzpwB/KnsXjJWb+WzCcLO6R/kqNjAilbSyMCMgTSdBoSGDiSi795GNLhvhjPT+cO+d3mGJ//Xn1n1U3VV3D4BmZ+1JDJTmKWrI4EV/Q+4HbhEOgWAZWUkB2uKGSHnwD86qrbqv5wANsAf0Akpl6Z+2Rp8YEdRBpAy9WWbfLZ4O1PsqNk5iXnRbfuTv604Lmhpb1gw7gjRXtqhiLGZyGews4dhfiYZx4hTSPUa6CPYTEays6+GTX6tBn8jq8f476+dr54LcVhDLw8FghFDp+v0VKXL5tVxwyEAnYLl80Fa1FLO/nWwSu0FfxBXj9W9KoKMZUVNVBxez5oQXipXEs1YDCxxIGMzawmtfeFkz8K+IO9fixLoZklswyQGyO1XflADnp+OqHDuQBmkYk4+BHRBu9HXQVnfxTwd7D4gaB2oxOjvOq4c/kAboC65hgY7KVQ2kGNPo2C5NGhdkJc+pU4ekVK5A8S43d8/djypYrYf8UQvJIw0eldAp+OpbGc2TbM5ZkNB9iDf8aYMnn17FdckV0z+83gV1fdI/0LdOUrCGYr8/cYE2BEgOBTHhfxVizYNrF1Ni4g46MNsQEluqZubZ/+Xjoh+PIlJE9PKeBhlkV8OVKkboBmcjgqIo4NWiAEGc31/snh2uEA9lzAmwV+l+RuIV9Y78hF9yBZdNuwmVFnyR4fAizcwyCgfl3k7qmm9rl3eYATgB+vuksjP2yqmzidbIdk84HRofYNFkwNvgCeGj0dXV3BvM/170f+hOAbMV7rBonZiB87SNBGmp9BPo97A+SHkOJFVl1yIa2eOfi6tp85gsfYL8ySBOsIUTtS5wDHMVeAm2jEqI0+uiloB0rtRt4A+ByRuOgQ48nuSMGT5C7ZcyA1I251iR7Hm+iYsi3TbvndVvD7vX7sprfzrhhGXzSlFpm2UNOk4NHTwySAJHMRCmGiCBXGGGG80OZjCghNZ6+dzlHBz379WHbusIhTaKE/3MPRuJbOrHfhg+kgSd1iAk/+Ik0HrAUGCt+XbVZ6r5/OTuS3u/qprx/Lz10FT82uWRHz9bd9yYuThAzTebB0FrpxEIB751kCDBVyMPQER0t+FpWThwswXweL8cq5p0Uc5agO8BiumT0CG1GKOHhKMeOaN069fPLuETKArohO+ywuwHydAXzFaefkAjskIHq6NSS/LcoiXQo/8vkdj/up1Sb4rK/K2atDZN53yVg99d13z+rz81IGfNF6BIwSeLUS9AnCuClvOhRiGjwLZ3YAgkc30wktO7td3l1zNPCKNPKFK6jnS2hsbH+KwdR1FyyefI4R0kCOGB1CTIODxnitd/qpOPhFAzEuvyy5Oxd/goNGfJAEEPNF7pgV4gCAvEHJ3rFP5BN0QwOZjQVyBWyB7yRPjbIrnhLrTsUhGJMUP8KIuH3kvgJpQ0LP0odAS6L9U39A+6RfAd7vGToB+F6T10wrO4Y5fcjhAWsiR0J2xO0JFQ4a7sRJpFBxpoNFehh4LeuvwIh6wL9+DO++qndet1Q9oL7z7gAvDuFTc2TAxgAciIdEHC0cmQqeNJ+Dj+Two4G/PFZ3WZ97RvAlcpGPCDpPI5NBavS0uqB3APlCo9xhlLq9B/kO8Dfg3z48DnwveeFT9cwukD+SWxfzODL9YtGAjSscFxg+AhkVeQdIOe5FWPDfcAFG1Gnxb3r5668HB691INA/MWbDNjp1g4QuEgMPmAuwcQU+HdN47sO1Dkhn0CaeNJ18T4z/8el6H+ZF+bbMhqqH1D5v5drJ644M+UHgnwvgaRzmLRIPgVWFjCvtY6BdrTh3oYeA36fqITWndNolxBgKCy3wN4bluGBO4ImBy3SbQOURnjh92hFRZtmGC4k9po5XoPPATnWB//G/l//+/MfjLL75DJZqO9J4b1QJyIhxXIvxzMPDILk1B8s+ge3gHcGBxeb1a67NQ8Bf77m/vCvcfl1ZdW8tyyVrLNnrThPMmXDEg2l+pf5TzTyNFQAfmZln7j2WAvuYHgL+jfx/fiw9b7G26v5KrhUVVuy5Ieml+ZVGtxzY3TTln7BmSPGZxQtXQAIAdoFnleKITk0m3xvjvw1znwZ+uUTwfmmYPYXFlmrFg3LhmX9PPVVjOYKnAQFqhPyBFojEjbAD5Xk1LpAcGfcHX3msbkvVvZUks1zyOxZnJfeiGWWRlidZlBtvgmYGgVEk8Z94iAjTu7QODMsESoeqJy4PmEu+1+LnV91bC+RCOfhOMwLbpLsRIm8qAm505IGShCUfyhicSAoIdIDo424z+OFgkdVf0ePBkysgvzybG7zaKqJhRXFEpQ1pN7F68PQxgQ+k8Sw5hFHCwdPBQxrsPO0kAb4dLJr1V3QA8HiJml+XlyZNtvIrH1KQx5I0EAfCGzw1ifNR2DU69WU0BaHIA//YaYPEM1hj1WXV13cfATyqGeUYyjwkhJROowlGuoMOm4g4WaQGkPSvyJ16gDaoBwm0f2v8tIPvbV0LCcnp0wDPEvYIhfD2K4R0OA6rYaMq0jGRekAiQWEwrjr9NWOHlK7vPhb4bU9hke0iwvNbbZQasXZq8clxRF4NwmbgedbIPzXslsOdmdc/LfiSQRTAg5ViJCB0yEBIU3VoIzDgpAWexbMKe8GLvXbBj5862DHdwlO5SLEEOijk2OFRHVIGjOFQeVaS1Y1VU++vDFQ5LCaSPw/464UZPfPblZPXDxCx2pewLbyB7AF+gghPwfOP1MmngEFrDuyI9CdpR3TcIvjlGgyeegW8fiwN4zgLEKVIrWwqT4K/CPNl8HSA8DsHZDuf0E7R+cAPnbsOHqCwLSQI0A83o5bLAFH+ndrjGYJw7SxnxJPKnIfoLW6YR/5s4JvkBdEsxqMVY6Smpon2jbYP+X46hvkC1iLG+Iwhb5kNmdx5aKceTYKH61snr9y5ySqBCfpSQPMJfNIeUs2BoIUWufdeinB3L5PLIloe47MWLIIH1cEHtEitS8uCDAEfSuCBOrsFHxJOFs95QerxWa4XVZsGB1E+Kdg7jfxm8K8fw4eXfX62vKTayaMLRk9Ou0TgpeMhvEcW/JllB8SZ6oQcLmbgMZLkHh+iBuzOB504I7rzMOB//vE1vnx44/+3HX7EuKCqsyeuOXUi8H2BXkzhtjGQs0EApWEEJD8B8ZwmB9lcMQrw7HMdvNh7GPCXb1m8/Db9ly3rrq9NPmrgg1Ivv65kXQ4zNPIXRmAYBelPrAftPT8d7tVXgJ9G/pgW3/B+HU9b55Nv8Zl6dPjAXQJL4Aho5i2S8y60Uz2vlAJU09EowdeNolczYvxv02N8C/z4uytEGpY3kGJvSBNqdts2gWfTtbQ1ykhfzNYWE6eOROtsbVvtV8BGtBn86qpbRevFu02eVshMuAmejJhb3kCg45YsxkPzUFw2Ujy59mk/O/hGjI9N8HAJEz7RnxJ4TMOTScNfknhk3GX32Sji57UNvPYrYGu0Hfy36/l9Hq56o+rk0yWkF5rN7nNajAnyjzECQJ7hiTuyvK5IhodobjP4Y8T4b2/IXz7Hbx9Gq96o1gIeDcyRYq2nUjyfYzM2hjEEWn0GnoQF4jciKdtovqxJaf1W8Ndp3FtKv9OLCgq6XKCGycOqjJg9dVlVqoFaOXXcCVDJc5ORIlZwupsvqv9tfFVtBX+ZzsXvH/Z5C1VJt6vX8QVabl3gc/t6Bg6CjKHAuEc2JGhJYvBYFf/UMP6iJv342WZXf/mezS9/fiffqCymPIM9K1eQwHckeEquHdST1o5MDiJgQW7jxAvQE8YwEGTt1OnH2D8KUUcBv7rqjtLFUHj71yCvB9/CgBJROKKXj8TcMXYw++XevBL7Y5ZFOHitdNnkb1evRV6ZZ5VqZVZMJwIR/mJhuwq+MWvDJhvg1TF6kBifND256xrMAz99mH0qNIfWijk97qMZW8Qcgh4U6ZEY7JXmWwGxcAnmpPWHtfi+wTzwQ2iN6jBRFzyBHwkUfGfWX+bOSXKQJQMd3c43Pzv4LvU9gncz4r4LzmdgMABw9MAIIJ4hlV2S/kgTOIo9iwmtbh8Y/I7vj+9RF/nkckn6VjoSvC9wxpw+0tKQ6kVBWW4ibkPsbY5E/YBjgN/1bdI96iMP+K5/xVLn6B5M4QX4tBscAUSJmKd1NMbLjKHYjbqOAX7X98c3FUJ7Og9HquCTfcP/2IJPIMqrA1dAn93AAIENEO8RMcifG/xDLf526Ua8/e0zybKXXZH/I8eFqGQHzNejkRNzJnVESBtDpTE8vKVjgFfeH99Z9QQtl66HvHTiBRbc4lN4zv18xGERCHCI+pgJkECwuAb09RnlLidwEPCrq56gdJ26yJPkO0KUpSzowg2aM6bvWcN4BNZG4jnL/TXw5ROqysEjzbEvVtH0isV4ZrsJKbvbAs3SA6UvJ0MpsAZ4QpGfhoMf1iB5moanbfwzJOPUfvmhJCEUdfDgEfGjHsJJ7SeK8aurnitBXubVBWn+PxJDDosXUMHrpkujBZZa/ElUrwpuLo4NKgfP9D67azUyaxJXHPJ5jAnEsoP4l1VF99NOtMCnAVfvq4PnouQlmlZXSsewGSAeow+pzMvQQjT+681ANungB0XIrwKvOdk8FtCxwGrmsUJMF/Lq1C44+FWi5KWxaVKzMn5EpDCptw4w4+PH8il9HWVe3GP8OvEMT72EpXRO2m9+iJitZcv3NAtk6ZpSb155XypylYPPNfQGExmiYwM83aI9u4fpIEzza4G9liJU5eAVdT2AmX0mmOixZL2Vb6Vw87p5Vq8R5Yn/qBy8ppEHMCnvPBsghyrGqARjyBD1tJJkFOgQus2cyMGrGn0Ak6RfpaDba5pyIPFRhn8oNY+MAAev6/3wwv0s8HhwtpjbAK/dCCpUH4KDL2iMPAUV9WwLgnqjJhG5i1k7evhh8JdjHHxRY+Qr6ytoum2zh3xfA8/DOQv3t/91ensHX9fg71+WRHP9IP7WD6YHqM691FJlvi8bcfAVzSFP4n9u0PLYfIUuXxlsNNZxyTzGNzSFPAWfh3D+bEbVVfeDv9eK7dOCH03udeUOGCdgkdFupGeNG0Dwd8cIcfANzSIvVvVYzEfcFfAa4FjwAiWHQlR+4qQqbZ2xohODn+Hug/wnbvEQ3Mlei3WU6tQPTilFzpXfiNIq6enFU4PfTl5AYu79tpSPdl6yvSHwYhKhHLPyHWS2wG8mr4FfYONyOzlcvREcS1lCzUs7+G0qBvrO2BiII0/WHtif9OConnQhxve03ATvMb6kAvluSxHFwL2DuVODF5WuufWmNMh0mF+vXl/1vaSSXws+pkT/9oln8qLSehBfKQffL83oW4G2tDnirbzk9IulGHi58rtWDn5AKvnKZLoyLtiKDv8xHJnpsbqxQti4yvRNv5NmhYrZfQO8spvhZPdjag4gB1+c91fl4AdVyu7Hweurr1EDr7bDbsR29Z3KwY+qSL4a43sjchO8FjNW+PojvWny9ePnn18uv2o6WvW9pWf3jYvfu1+L8dXCvT/CxXUg8JefQrn8cvk930K1UorRr8uxW/fZK+6CZoPnBv/j7/+8/oD1Y378aEw5+W5XXpiuq+WrmeGmad2h3h//Zu7ff7v+cvlg1Y+QRN95+flhW8BjGvBgg5+Q3N1eTZJzPyL4nHz3s60q+HKKr+9cmcsvmsndTFYP6nk+Q7uhQWJyCPXUrPYMlj577xt/Dn6T2uRzH40PwcI0fOWae9edW11TuU9z9Xd/C9UGtdBnHPBpOPpcnB6rgxbDG6PklOAf9RaqLWqQlxxunj7SNRqM49qxcnsLbBf4udxnTOce8BaqzaqjlwsxsNCG4JOrnwK+K8YfDPxD3kI1QSPP4NIXlEhPPgd8V4c3ViC02dXv9haqvTWAHh6DYjfcgtxCNw/F+A5N+lIYyl5WDxohH/O1mHzTfpr1vSlaZX33M4PvmtkV19nuCf7y0+yz25oF/mTJ3aL3g/l9367perP3w4Ifr/oYqqOv0b1jEnPpooOfrfcV9m2zvgf+95PeHM60GfyD30I1Q0X07bnJsMN/9E057Eh997HfSTNPKno+F1eJjYKnx/eOgV24z1m5i2dN7ogaz+cUCG8A31t0+gwem6/IisVflKHvAD/qusfBb/zyX7UnFR37LVSTlZNXHozvlp4gtJ7Xa3VpfX+UCioykNVTVeZ2bQuT36NRXk/JDuyw2S3zjYYcPFNrRacigaIEnn0Jq92bztbG5eCFVqOHe3Lpfq0K/nK3J3Rem+p3QB4e41dXfVjVVnQqWsDjnXo9xjPwNXrYifJ3MTfIwStahx6deOXMmbHWjiXc+VEOfletNHsgWzx5uqvCUHmzUrvQiBx8Se9X2313+lY8jjUsj5pze8DBlxWG2JNb9xXwIl6rDGWju9wIcvBFXfF1o+9Z5uvy0xtmlCNy8EUlSH0kGNKijTbB3wm7g68IIfWYfV9cbxy1IrKslYMvi1zbDvSd379szNw3PP45Jgffq5VZ/lj93Uc7+DvqPWinugeOd/B31jT4y/fwSBKxovyW9uu7HbymCfRvC3xxdLFgnhz8Sm2Ev4DfL3b0dKAiB1/VevohPhD60oGaHHxT6+CnMg+7SA5+it5nah5wg+7gzy+NbVlx2g3WdXLw+6jlAS5y8Fb1yB+PcPBG5eCNysE/WI9y9w7+sXpYgufgHysHb1QO3qpGv2c9KyVw8KfSPAfh4E+lo4F//TX/XQQHv4OOA375QQztJzEc/A46Toz/8ekNuVv86TTB1f/49Mu/HPzZNCXGv35UfvvIwR9antUblYM3qu3gz/cWKlecAP6Mb6FyTZjHn/MtVK6t4M/6FirzmrCAc9K3UBmXZ/VG5eCNahZ4T+5OJrd4o3LwRrUZ/BO8hcqkNs/jn+edNLY0ZeUuenJ3OrnFG9WUlTuP8eeTZ/VG5eCNysEbGBsyQwAAA2ZJREFUlYM3KgdvVA7eqBy8UTl4o3LwRuXgjcrBG5WDNyoHb1QO3qgcvFE5eKNy8Ebl4I3KwRuVgzcqB38APeIr5Q7+8XrIT5c7+MfLwRuVg7cqj/Guu8nBG5WDNyoHb1QO3qgcvFE5eKNy8Ebl4I3KwRuVgzeqPcG7jqz9wMMAOFPJU3V2x5IO3mhJB2+0pIM3WtLBGy3p4I2WdPBGS84A7zqhHLxROXijcvBG5eCNysEblYM3KgdvVA7eqBy8UW0A//rr7Y1krx+v7yZ7Ceo7yjpK/vgU8pfaTW3z8prU2+vzlra6m8xKdp8mlly6vaLNpeSKNpdP5TbXg/++9OXH71/jy1vt3z63SmQlP32+lLy80vDlw5qSvW1euvj671/j7fWJb211N5mV7D5NLLl0e0Wb6YTH21w+VdpcDf7bu/+7Wd/r3/68vH/05x9fV5X8/evl5aXJksdK9rb5/XLu16u3tNXdZFay+zSx5NLtFW0uJVe0uXyqtLnd1S8Wf30hZefIfGX4lv+tKDnS5q2FFU3ykiNNYhOXbq9ocym5rs3GeU6I8bc4cnExvUMzDZmLw3739fLW2kHwS8mBNn9++e3yv6WtkSZ5yZHTXEou3V7R5lJyXZuXT5U2t4O/dOv7knp0BiOSov33H2ssfinZ3+aPT7frMW7xomR/k1AyDlu8KLmqzeunXS2ejqpB8PEaOvtjvCzZ3ebrx89YaCTGZyW7m8SSS7dXtBlHwWPJ26ddY/xi8Rf+P/8xhO+aHXy4+qTOrF6U7G0Tr8fSVneTWcnu06T0rt1e0WbEIDHY5vKp0uZG8Nd+hfDu63Wy+a4zdSElhybVWcnONl+u3yz5vCRKI03mJXtPk5QcnMfnJcfbXD7tMo93nVoO3qgcvFE5eKNy8Ebl4I3KwRuVgzcqB29UDt6oHLxROXijcvBG5eCNysEblYM3KgdvVA7eqBy8UTl4o3LwEb+s0PuQ9zPIwSP4773f9n0GOfgr+J9ffvkTvpJpQg7+Cv7b9ZtH7upt6ecf/7F8P87Bm9LPL//1P9evmzh4W3pz9S/u6g3qktxdE3sHb0sX6lfmDt719HLwRuXgjcrBG5WDNyoHb1QO3qgcvFE5eKNy8Ebl4I3KwRuVgzcqB29UDt6oHLxROXijcvBG9f86tMJHNbYKiAAAAABJRU5ErkJggg==" alt="plot of chunk unnamed-chunk-1"/> </p>
+
+<pre><code class="r">
+# percentage of parameters from simulation within the 95% ellipse
+tmp=rep(0,length.out=nsim)
+for(i in 1:nsim){
+  tmp[i]=(parsim[i,]-as.numeric(parest))%*%Sinv%*%(parsim[i,]-as.numeric(parest))
+}
+sum(tmp &lt;= qf(0.95,2,58)*2)/nsim
+</code></pre>
+
+<pre><code>[1] 0.957
+</code></pre>
+
+<pre><code class="r">
+# session Info
+sessionInfo()
+</code></pre>
+
+<pre><code>R version 3.0.1 (2013-05-16)
+Platform: i386-w64-mingw32/i386 (32-bit)
+
+locale:
+[1] LC_COLLATE=English_United States.1252 
+[2] LC_CTYPE=English_United States.1252   
+[3] LC_MONETARY=English_United States.1252
+[4] LC_NUMERIC=C                          
+[5] LC_TIME=English_United States.1252    
+
+attached base packages:
+[1] stats     graphics  grDevices utils     datasets  methods   base     
+
+other attached packages:
+[1] minpack.lm_1.1-7 deSolve_1.10-6   reshape2_1.2.2   ggplot2_0.9.3.1 
+[5] knitr_1.2       
+
+loaded via a namespace (and not attached):
+ [1] colorspace_1.2-2   dichromat_2.0-0    digest_0.6.3      
+ [4] evaluate_0.4.3     formatR_0.8        grid_3.0.1        
+ [7] gtable_0.1.2       labeling_0.2       MASS_7.3-26       
+[10] munsell_0.4        plyr_1.8           proto_0.3-10      
+[13] RColorBrewer_1.0-5 scales_0.2.3       stringr_0.6.2     
+[16] tools_3.0.1       
+</code></pre>
+
+</body>
+
+</html>
+
Add a comment to this file

04_learnR_parfitODE/concProfileData.jpeg

Added
New image
Add a comment to this file

04_learnR_parfitODE/confRegion_simul95pct.jpeg

Added
New image

README

-This is a place to put codes and material I am using to learn different aspects of math, computation and visualization. I put some of these in the blog http://notesofdabbler.wordpress.com
+This is a place to put codes and material I am using to learn different aspects of math, computation and visualization. I put some of these in the blog http://notesofdabbler.wordpress.com

R_codes/20120804_CompApp.Rmd

-# Learning R using a Chemical Reaction Engineering Book
-
-[Chemical Reactor Analysis and Design Fundamentals](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/) by J.B. Rawlings and J. G. Ekerdt is a textbook for studying Chemical Reaction Engineering. The popular open source package [Octave](http://www.gnu.org/software/octave/) has its origins to the reaction engineering course offered by Prof. Rawlings. This book is accompanied by Octave and Matlab code for solving typical problems encountered in Reaction Engineering. 
-
-I figured that maybe one way to learn R is so see whether I can code some of th examples from this book in R. I am by no means suggesting that R can replace MATLAB/Octave for engineering problems but merely it is a way for me to learn the language.
-
-I started with the [computational appendix](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/web-appendix.pdf) listed in the book's website and am trying to work through some of the examples there. It will be good to refer to the computational appendix to follow the R code below.
-
-### Setting up a stoichiometric matrix, reaction rate vector and determining the rank
-
-```{r}
-# stoichiometric matrix
-stoi=matrix(c(0,1,0,-1,-1,1,
-       -1,1,1,-1,0,0,
-       1,0,-1,0,-1,1),
-       ncol=6,byrow=T)
-stoi
-
-# rank of the stoichiometrix matrix
-rank=qr(stoi)$rank
-rank
-
-# reaction rate vector
-r=c(1,2,3)
-```
-Given the reaction rate r=(r1,r2)', the rate of change of species concentration R is given by
-$$latex R=\nu^Tr $$
-where $\nu$ is the stoichiometrix matrix
-```{r}
-# rate of change of components
-R=t(stoi)%*%r
-R
-```
-### Example A1: Estimating reaction rates
-
-The stoichometrix matrix $\nu$ is input below
-```{r}
-# stoichiometry
-stoi=matrix(c(0,1,0,-1,-1,1,
-              -1,1,1,-1,0,0),nrow=2,byrow=T)
-stoi
-# number of species and number of reactions
-nspec=ncol(stoi)
-nr=nrow(stoi)
-
-nspec
-nr
-
-# true rxn rates
-r=c(1,2)
-r
-
-# true component rates
-R=t(stoi)%*%r
-R
-```
-
-### simulate 2000 measured component rates
-
-Add random noise (normally distributed with mean 0 and standard deviation 0.05) to true species rate vector R
-$$ R^m = R + \epsilon, \epsilon \sim N(0,0.0025) $$
-
-```{r}
-## simulate 2000 noise estimates
-e=matrix(0.05*rnorm(2000*nspec,0,1),nrow=2000,byrow=T)
-Rmeas=matrix(rep(R,2000),ncol=nspec,byrow=T)+e
-```
-The least squares estimate of reaction rate vector $\cap{r}$ is
-$$ \hat{r}=(\nu\nu^T)^{-1}{\nu}R $$
-
-```{r}
-## estimate reaction rates
-rest=solve(stoi%*%t(stoi),stoi%*%t(Rmeas))
-```
-I was trying different plot features in R and applying to this data of estimated rates
-
-I found the following function that plots scatterplot with marginal histograms
-```{r}
-# plotting scatterplot with histogram
-# downloaded from web
-#  http://www.r-bloggers.com/example-8-41-scatterplot-with-marginal-histograms/
-#
-scatterhist = function(x, y, xlab="", ylab=""){
-  par.default <- par(no.readonly=TRUE)
-  zones=matrix(c(2,0,1,3), ncol=2, byrow=TRUE)
-  layout(zones, widths=c(4/5,1/5), heights=c(1/5,4/5))
-  xhist = hist(x, plot=FALSE)
-  yhist = hist(y, plot=FALSE)
-  top = max(c(xhist$counts, yhist$counts))
-  par(mar=c(3,3,1,1))
-  plot(x,y)
-  par(mar=c(0,3,1,1))
-  barplot(xhist$counts, axes=FALSE, ylim=c(0, top), space=0)
-  par(mar=c(3,0,1,1))
-  barplot(yhist$counts, axes=FALSE, xlim=c(0, top), space=0, horiz=TRUE)
-  par(oma=c(3,3,0,0))
-  mtext(xlab, side=1, line=1, outer=TRUE, adj=0, 
-        at=0.9 * (mean(x) - min(x))/(max(x)-min(x)))
-  mtext(ylab, side=2, line=1, outer=TRUE, adj=0, 
-        at=(0.9 * (mean(y) - min(y))/(max(y) - min(y))))
-  par=par(par.default)
-}
-
-# scatter plot of reaction rates with marginal histograms
-scatterhist(t(rest)[,1],t(rest)[,2],xlab="r_1",ylab="r_2")
-```
-There is library cars that has a command 'scatterplot' for plotting scatterplot with box plots and has several other options
-In the plot below 50% and 90% ellipses are overlaid on the data
-
-```{r}
-# scatter plot of reaction rates with marginal box plots
-library(car)
-scatterplot(t(rest)[,1],t(rest)[,2],reg.line=FALSE,smooth=FALSE,ellipse=TRUE)
-```
-I tried ggplot2 also for the same plot
-
-```{r}
-# 2d contours/density with ggplot2 for reaction rates
-
-# create dataframe of reaction rates
-rest_df=data.frame(t(rest))
-names(rest_df)=c("r1","r2")
-
-library(ggplot2)
-ggplot(data=rest_df,aes(x=r1,y=r2))+geom_point()
-ggplot(data=rest_df,aes(x=r1,y=r2))+stat_density2d()
-ggplot(data=rest_df,aes(x=r1,y=r2))+stat_density2d(aes(fill=..level..),geom="polygon")
-```
-
-

R_codes/20120811_CompApp_pt_2_3.Rmd

-# Learning R using a Chemical Reaction Engineering Book (Part 2)
-
-In case you missed part 1, you can view it [here](http://rdabbler.blogspot.com/2012/08/learning-r-using-chemical-reaction.html). In this part, I tried to recreate the examples in sections A.2.2 and A.2.3 of the [computational appendix](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/web-appendix.pdf) in the [reaction engineering book](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/) (by Rawlings and Ekerdt). 
-
-## Solving a nonlinear system of equations
-
-This example involves determining reaction equilibrium conditions by solving the following system of nonlinear equations.
-
-$$latex
- \begin{aligned}
-   PK_1y_Iy_B-y_{P1}=0, \\
-   PK_2y_Iy_B-y_{P2}=0
- \end{aligned}   
-$$
-
-The relation between the variables $latex y_I,y_B,y_{P1},y_{P2}$ and extent of reactions $latex x_1,x_2$ are:
-$$latex
-  \begin{aligned}
-    y_I=\frac{y_{I0}-x_1-x_2}{1-x_1-x_2}, \\
-    y_B=\frac{y_{B0}-x_1-x_2}{1-x_1-x_2}, \\ 
-    y_{P1}=\frac{y_{p10}+x_1}{1-x_1-x_2}, \\
-    y_{P2}=\frac{y_{p20}+x_2}{1-x_1-x_2}
-  \end{aligned}
-$$  
-
-Here I have used R package [rootSolve](http://cran.r-project.org/web/packages/rootSolve/index.html) for solving the above set of equations. The library is loaded and the functions to be solved are defined in the R function fns.
-
-```{r}
-# load library rootSolve
-library(rootSolve)
-
-# function defining F(x)=0
-fns=function(x){
-  K1=108; K2=284; P=2.5
-  yI0=0.5; yB0=0.5; yP10=0; yP20=0;
-  d=1-x[1]-x[2]
-  yI=(yI0-x[1]-x[2])/d
-  yB=(yB0-x[1]-x[2])/d
-  yP1=(yP10+x[1])/d
-  yP2=(yP20+x[2])/d
-  F1=P*K1*yI*yB-yP1
-  F2=P*K2*yI*yB-yP2
-  c(F1=F1,F2=F2)
-}
-```
-Next, an initial guess of (0.2,0.2) is set for the variables and the equations are solved using the function multiroot (from package rootSolve)
-
-```{r}
-# initial guess for x
-xinit=c(0.2,0.2)
-
-# solve the equations
-xans=multiroot(f=fns,start=xinit)
-
-# object returned by multiroot
-xans
-
-# solution to the equations
-xans$root
-```
-The solution to the equations is accessed from the variable xans$root which in this case is (0.1334,0.3507)
-
-## Function Minimization
-
-Here the function to be minimized is
-
-
-$$latex
-  G=-(x_1lnK_1+x_2lnK_2)+(1-x_1-x_2)lnP+(y_{I0}-x_1-x_2)ln(y_I)+(y_{B0}-x_1-x_2)ln(y_B)+(y_{P10}+x_1)ln(y_{p1})+(y_{P20}+x_2)ln(y_{P2})
-$$
-
-The following constraints need to be satisfied for $latex x_1$ and $latex x_2$
-
-$$latex
-  0 \leq x_1 \leq 0.5, 0 \leq x_2 \leq 0.5, x_1+x_2 \leq 0.5
-$$latex
-
-First I used [constrOptim](http://stat.ethz.ch/R-manual/R-patched/library/stats/html/constrOptim.html) function for minimization. We need to specify the function to be minimized
-
-```{r}
-# function to be minimized
-eval_f0=function(x){
-  dg1=-3.72e3; dg2=-4.49e3; T=400; R=1.987; P=2.5
-  K1=exp(-dg1/(R*T)); K2=exp(-dg2/(R*T))
-  
-  yI0=0.5; yB0=0.5; yP10=0; yP20=0;
-  d=1-x[1]-x[2]
-  yI=(yI0-x[1]-x[2])/d
-  yB=(yB0-x[1]-x[2])/d
-  yP1=(yP10+x[1])/d
-  yP2=(yP20+x[2])/d
-  
-  f=-(x[1]*log(K1)+x[2]*log(K2))+(1-x[1]-x[2])*log(P)+yI*d*log(yI)+
-       yB*d*log(yB)+yP1*d*log(yP1)+yP2*d*log(yP2)
-  
-  return(f)
-}
-```
-The constraints need to be specified in the form $latex Ax-b \geq 0$
-
-```{r}
-#  constraint
-A=matrix(c(-1,-1,1,0,-1,0,0,1,0,-1),ncol=2,byrow=TRUE)
-A
-b=c(-0.5,0,-0.5,0,-0.5)
-b
-```
-Next, the function is minimized using constrOptim (starting from an initial guess of (0.2,0.2)). Here Nelder-Mead method is used since BFGS method requires specifying the gradient of the function by the user. R taskview [optimization](http://cran.r-project.org/web/views/Optimization.html) lists other options.
-
-```{r}
-# initial guess
-xinit=c(0.2,0.2)
-
-# minimize function subject to bounds and constraints
-xans2=constrOptim(theta=xinit, f=eval_f0, grad=NULL, ui=A, ci=b,
-            method = "Nelder-Mead")
-xans2
-```
-The solution can be accessed from xans2$par and is (0.1331,0.3509)
-
-Next, I also tried function minimization with ConstrOptim.nl from the package [alabama](http://cran.r-project.org/web/packages/alabama/index.html). Here the constraints are specified in terms of $latex h(x) \geq 0$. Even if gradient is not supplied, this function will estimate it using finite-differences.
-
-Definition of constraints in the format for constrOptim.nl
-
-```{r}
-# load library alabama
-library(alabama)
-library(numDeriv)
-
-h_ineq=function(x){
-   h=rep(NA,1)
-   h[1]=-x[1]-x[2]+0.5
-   h[2]=x[1]
-   h[3]=-x[1]+0.5
-   h[4]=x[2]
-   h[5]=-x[2]+0.5
-   return(h)
-}
-
-xans3=constrOptim.nl(par=xinit,fn=eval_f0,hin=h_ineq)
-
-xans3
-```
-The solution can be accessed from xans3$par and is (0.1332,0.3508)
-
-Since this is a 2 dimensional problem, the solution can also be visualized using a contour plot of the function.
-
-```{r message=FALSE,warning=FALSE}
-# Region of interest: 0.01<=x1<=0.49, 0.01<=x2<=0.49
-x1=seq(0.01,0.49,by=0.01)
-x2=seq(0.01,0.49,by=0.01)
-
-# vectorizing function eval_f0 so that it can be evaluted in the outer function
-fcont=function(x,y) eval_f0(c(x,y))
-fcontv=Vectorize(fcont,SIMPLIFY=FALSE)
-z=outer(x1,x2,fcontv)
-
-# filled.coutour and contour are overlaid with the minimum point (0.133,0.351)
-filled.contour(x1,x2,z,xlab="x1",ylab="x2",
-               plot.axes={axis(1); axis(2); contour(x1,x2,z,levels=c(-3,-2.5,-2,-1.5,-1,0),vfont=c("sans serif","bold"),labcex=1,lwd=2,add=TRUE); 
-                          points(0.133,0.351,pch=15,col="blue")})
-```
-
-MATLAB/Octave functions for solving nonlinear equations (fsolve and fmincon) have been used in Chemical Engineering computations for a long time and are robust. R has traditionally not been used in this domain. So it is hard to say how the functions I have used in this blog will perform across the range of problems encountered in Reaction Engineering.
-
-This post was created with R 2.15.1 and knitr 0.7.1 (converted markdown to html)

R_codes/20121209_CompApp_pt4.Rmd

-# Learning R using a Chemical Reaction Engineering Book (Part 3)
-
-In case you missed the previous parts, the links to them are listed below:
-* [Part1](http://rdabbler.blogspot.com/2012/08/learning-r-using-chemical-reaction.html). 
-* [Part2](http://rdabbler.blogspot.com/2012/08/learning-r-using-chemical-reaction_11.html)
-
-In this part, I tried to recreate the examples in sections A.3.1 of the [computational appendix](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/web-appendix.pdf) in the [reaction engineering book](http://jbrwww.che.wisc.edu/home/jbraw/chemreacfun/) (by Rawlings and Ekerdt). 
-
-## Solving a system of ordinary differential equations
-
-This example involves reaction (Benzene pyrolysis) in a plug flow reactor. The actual reactions happening are:
-
-$$ (Rxn1) \; 2B = D + H $$
-$$ (Rxn2) \; B + D = T+H $$
-
-
-The rate of each reaction is given by:
-$$ r_1=k_1(c_B^2-\frac{c_Dc_H}{K_1}) $$
-$$ r_2=k_2(c_Bc_D-\frac{c_Tc_H}{K_2}) $$
-
-The feed to the reactor consists of 60kmol/hr of Benzene (B). The temperature of the reactor is $T=1033K$ and the pressure is $P=1atm$.The rate constants and equilibrium constants for this example are:
-$$ k_1=7\times 10^5\; L/mol.hr,\; k_2=4\times 10^5 \; L/mol.hr,\; K_1=0.31,\; K_2=0.48 $$
-
-```{r, tidy=FALSE}
-# set working directory
-setwd('~/R/RxnEngg_Rawlings')
-
-# load libraries
-library(deSolve)
-
-
-# Appendix A.3.1: Solution of Differential Equations
-
-# Benzene pyrolysis example
-
-# Parameters
-# NBf - feed benzene rate - mol/h
-# R - Universal gas constant
-# T - Reactor temperature K
-# P - Reactor pressure atm
-# k1 - rxn1 forward rate constant L/mol.h
-# k2 - rxn2 forward rate constant L/mol.h
-# Keq1 - rxn1 equilibrium constant
-# Keq2 - rxn2 equilibrium constant
-pars=c(
-NBf=60e3,  
-R=0.08205,  
-T=1033, 
-P=1, 
-k1=7e5, 
-k2=4e5, 
-Keq1=0.31, 
-Keq2=0.48  
-)
-```
-
-The governing equations for conversion versus volume in a plug flow reactor is based on extent of each of the reactions:
-$$ \frac{d\epsilon_1}{dV}=r_1,\; \frac{d\epsilon_2}{dV}=r_2 $$
-The initial conditions (corresponding to feed conditions $N_B(0)=60kmol/h,\;N_D(0)=N_H(0)=N_T(0)=0$) are that the extent of reaction is zero. 
-$$ \epsilon_1(0)=0, \; \epsilon_2(0)=0 $$
-
-The flow rates of each component along the reactor volume can be calculated from reaction extent
-$$ N_B=N_B(0)-2\epsilon_1-\epsilon_2, \; N_D=\epsilon_1-\epsilon_2, \; N_H=\epsilon_1+\epsilon_2, \; N_T=\epsilon_2 $$ 
-
-These are setup in a function that can be passed to an ODE solver. In this case the ODE solver we use lsode from the R package [deSolve](http://desolve.r-forge.r-project.org/). The inputs to the function are:
-* Variable over which the integration is done (Volume in this case)
-* The state variables of the system (Extent of the two reactions)
-* Parameters that are needed for description of the system (Rate constants, Temperature, Pressure, etc.)
-The output from this function is the rate of change as described by the equations previously.
-
-```{r, tidy=FALSE}
-# function that will be passed to odesolver
-# vol is the variable over which the system is integrated (equivalent of time in batch reactions)
-# ext is the extent of reactions 1 and 2
-# params are the parameters passed to the system
-rxnrate=function(vol,ext,params) {
-     with(as.list(c(ext,params)),{
-        NB=NBf-2*ext1-ext2
-        ND=ext1-ext2
-        NH=ext1+ext2
-        NT=ext2
-        Q=NBf*R*T/P
-        cB=NB/Q
-        cD=ND/Q
-        cT=NT/Q
-        cH=NH/Q
-        dext1=k1*(cB*cB-cD*cH/Keq1)
-        dext2=k2*(cB*cD-cT*cH/Keq2)
-        return(list(c(dext1=dext1,dext2=dext2)))
-     })
-}
-```
-
-Since the reaction start only after the feed enters the reactor, the extent of reaction is zero for both reactions at the beginning of the reactor (V=0L). The set of volumes where the concentration and reaction extent is computed is chosen in this case to be from 0L to 1600L at every 50L. The ODE solver lsode from [deSolve](http://desolve.r-forge.r-project.org/) package is used to solve the system of equations.
-
-```{r, tidy=FALSE}
-# initial extent of reaction (zero in this case for both reactions)
-extinit=c(ext1=0,ext2=0)
-# Volumes where the concentration is reported (in this case 0 to 1600L at every 50L)
-vols=seq(0,1600,length=50)
-# Solution of the set of differential equations using lsode solver in deSolve package
-extout=lsode(times=vols,y=extinit,func=rxnrate,parms=pars)
-```
-**extout** contains the extent of reaction vs volume data. That is used to compute mole fraction and conversion at different volumes along the reactor.
-
-```{r, tidy=FALSE}
-# Calcuation of mole fraction and conversion from extent of reaction at different volumes
-extoutdf=data.frame(extout)
-NBf=pars["NBf"]
-extoutdf$conv=(extoutdf$ext1*2+extoutdf$ext2)/NBf
-extoutdf$yB=(NBf-2*extoutdf$ext1-extoutdf$ext2)/NBf
-extoutdf$yD=(extoutdf$ext1-extoutdf$ext2)/NBf
-extoutdf$yT=(extoutdf$ext2)/NBf
-extoutdf$yH=(extoutdf$ext1+extoutdf$ext2)/NBf
-```
-Next conversion and mole fraction is plotted as a function of reaction volume
-
-```{r, tidy=FALSE, fig.width=10, fig.height=5}
-# load library ggplot2 for plotting
-library(ggplot2)
-# load library reshape2 for data reshaping
-library(reshape2)
-# plot of conversion vs volume
-ggplot(extoutdf,aes(x=time,y=conv))+geom_line()+
-  scale_x_continuous(breaks=seq(0,1600,by=200))+xlab('Volume (L)')+ylab('Conversion')+theme_bw(20)
-
-# plot of mole fraction vs volume
-tmp=melt(extoutdf[,c("time","yB","yD","yT","yH")],id.vars=c("time"),variable.name="moleFraction")
-ggplot(tmp,aes(x=time,y=value,color=moleFraction))+geom_line()+
-    scale_x_continuous(breaks=seq(0,1600,by=200))+xlab('Volume (L)')+ylab('moleFraction')+theme_bw(20)
-```
-Info on the session
-```{r, tidy=FALSE}
-print(sessionInfo(),locale=FALSE)
-```   
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.