himmelex

Himmelblau exercise 9.9

This is exercise 9.9 from David M. Himmelblau, Process Analysis by Statistical Methods, Wiley, 1970.

We model the reactions

A + B  (k1)-> C + F
A + C  (k2)-> D + F
A + D  (k3)-> E + F

The derivatives can be written as

dA/dt = -k1 AB - k2 AC - k3 AD
dB/dt = -k1 AB
dC/dt =  k1 AB - k2 AC
dD/dt =          k2 AC - k3 AD
dE/dt =                  k3 AD

The system is written in file himmelode.m and the sum of squares function in himmelss.m.

clear model data parama options

data.ydata = [
%  Time (min)     [A] (mole/liter)
            0     0.02090
         4.50     0.01540
         8.67     0.01422
        12.67     0.01335
        17.75     0.01232
        22.67     0.01181
        27.08     0.01139
        32.00     0.01092
        36.00     0.01054
        46.33     0.00978
        57.00     0.009157
        69.00     0.008594
        76.75     0.008395
        90.00     0.007891
       102.00     0.007510
       108.00     0.007370
       147.92     0.006646
       198.00     0.005883
       241.75     0.005322
       270.25     0.004960
       326.25     0.004518
       418.00     0.004075
       501.00     0.003715
    ];

Initial concentrations are saved in data to be used in sum of squares function.

A0 = 0.02090; B0 = A0/3; C0 = 0; D0 = 0; E0 = 0;
data.y0 = [A0;B0;C0;D0;E0];

Refine the first guess for the parameters with fminseacrh and calculate residual variance as an estimate of the model error variance.

k00 = [15,1.5,0.3]';
[k0,ss0] = fminsearch(@himmelss,k00,[],data)
mse = ss0/(length(data.ydata)-4);

params = {
    {'k1', k0(1), 0}
    {'k2', k0(2), 0}
    {'k3', k0(3), 0}
    };

model.ssfun = @himmelss;
model.sigma2 = mse;

options.nsimu = 1000;
options.updatesigma = 1;

[results,chain,s2chain] = mcmcrun(model,data,params,options);

Sampling these parameters:
name   start [min,max] N(mu,s^2)
k1: 14.4019 [0,Inf] N(0,Inf)
k2: 1.5663 [0,Inf] N(0,Inf)
k3: 0.290424 [0,Inf] N(0,Inf)

figure(1); clf
mcmcplot(chain,[],results,'chainpanel')
subplot(2,2,4)
mcmcplot(sqrt(s2chain),[],[],'dens',2)
title('error std')

Function chainstats lists some statistics, including the estimated Monte Carlo error of the estimates.

chainstats(chain,results)

MCMC statistics, nsimu = 1000

                 mean         std      MC_err         tau      geweke
---------------------------------------------------------------------
        k1     14.415     0.67092    0.053086       9.598     0.98736
        k2     1.5668    0.039761   0.0039085      11.468     0.99983
        k3    0.29106    0.012697   0.0012202      13.761     0.99665
---------------------------------------------------------------------

figure(2); clf
[t,y] = ode45(@himmelode,linspace(0,600),data.y0,[],mean(chain));
plot(data.ydata(:,1),data.ydata(:,2),'s',t,y,'-')
ylim([0,0.021])
legend({'Aobs','A','B','C','D','E'},'Location','best')
title('Data and fitted model')