Source

Blog_NotesOfDabbler / 04_learnR_parfitODE / compAppendix_parest.html

Full commit
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
<!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>