Bayesian Parameter Estimation (MCMC)

We are considering the boundary value problem

Where \(T_s(x)\) denotes the steady state temperature of an uninsulated rod with source heat flux \(\Phi\) at \(x=0\) and convective heat transfer coefficient \(h\). We fix the thermal conductivity at \(k=4.01\), the cross-sectional dimensions of the rod at \(a=b=0.95\) cm, the length of the rod at \(L=70\) cm, and the ambient temperature at \(T_{amb}=22.28\).

We are given experimental data on the state-state temperatures at various locations \(x\) along a copper rod, which we list in the table below:

Our goal is to use this data to build estimates of the parameters \(\Phi\) and \(h\) from a Bayesian framework. We have previously constructed frequentist estimates of \(\theta=\left[\Phi, h\right]^T\) as \(\widehat{\theta}_{OLS}=\left[-9.93, 0.00143\right]^T\). These values are the initial parameters we use in our Metropolis Algorithm.

Our jump distribution has a covariance matrix specified by \(V=R^TR\) where \(R\) is the Cholesky decomposition of \(V\). We construct \(V\) as \(V=\widehat{\sigma^2}(X^TX)^{-1}\) where \(X\) is the sensitivity matrix for the problem at hand, \(\widehat{\sigma^2}=\frac{1}{n-p}\sum\limits_{i=1}^{n}(y_i-\widehat{y_i})^2=\frac{1}{13}\sum\limits_{i=1}^{15}(y_i-\widehat{y_i})^2\), and \(y_i-\widehat{y_i}\) is the residual for the \(i^{\text{th}}\) observation; it is the difference between the observed data in the table above and the analytic solution given by:

Where we have abbreviated:

In MATLAB, We can create a target function to minimize the sum of squares like so:

%%%%%I. heatss.m%%%%%
function ss = heatss(params,data)

 udata = data.ydata;
 xdata = [10 14 18 22 26 30 34 38 42 46 50 54 58 62 66]';

 %%%1. Input Dimensions And Material Constants%%%
 a = 0.95;               %fixed, the cross sectional area%
 b = 0.95;              %fixed, the cross sectional area%
 L = 70.0;              %fixed, the length of the rod%
 k = 4.01;              %fixed, the thermal conductivity%
 u_amb = 22.28;   %fixed, the ambient room temp%
 Q = params(1);   %parameter of interest 1 (Phi)%
 h = params(2);    %parameter of interest 2%
 n = 15;                %number of observations%
 p = 2;                  %number of parameters%

 %%%2. Construct Constants And Solution%%%
 gamma = sqrt(2*(a+b)*h/(a*b*k));
 f1 = exp(gamma*L)*(h + k*gamma);
 f2 = exp(-gamma*L)*(h - k*gamma);
 f3 = f1/(f2 + f1);
 c1 = -Q*f3/(k*gamma);
 c2 = Q/(k*gamma) + c1;

 uvals_data = c1*exp(-gamma*xdata) + c2*exp(gamma*xdata) + u_amb;
 res = udata - uvals_data;
 ss = res'*res;

From the data in the table and the above formulas, our estimated variance is \(\widehat{\sigma^2} \approx 0.0412\) and our estimated covariance matrix \(V\) is:

We can implement the Metropolis Algorithm with 100,000 iterations of our chain with the following (the DRAM implementation is similar and omitted):

%%%%%III. heat_code_mcmc.m%%%%%
%Required functions: kde.m
%Required data: final_al_data.txt

load final_al_data.txt  
xdata = [10 14 18 22 26 30 34 38 42 46 50 54 58 62 66];
data = [66.04 60.04 54.81 50.42 46.74 ...                             %Copper Rod%
          43.66 40.76 38.49 36.42 34.77 ...
          33.18 32.36 31.56 30.91 30.56]; 
vals = [10:.1:70];                                                                  %A large range%



%%%1. Input dimensions and material constants%%%
a = 0.95;        %cross section of rod, fixed%
b = 0.95;        %cross section of rod, fixed%
L = 70.0;        %length of rod, fixed%
k = 4.01;        %new data given%
u_amb = 22.28;   %New data given%
h = 0.00143;     %Param 1 of interest, based on Proj 2%
Q = -9.93;       %Param 2 of interest, based on Proj 2 (Phi)%
n = 15;          %number of observations%
p = 2;           %number of parameters solving for%



%%%2. Construct Analytic Solution%%%
gamma = sqrt(2*(a+b)*h/(a*b*k));
gamma_h = (1/(2*h))*gamma;
f1 = exp(gamma*L)*(h + k*gamma);
f2 = exp(-gamma*L)*(h - k*gamma);
f3 = f1/(f2 + f1);
f1_h = exp(gamma*L)*(gamma_h*L*(h+k*gamma) + 1 + k*gamma_h);
f2_h = exp(-gamma*L)*(-gamma_h*L*(h-k*gamma) + 1 - k*gamma_h);
c1 = -Q*f3/(k*gamma);
c2 = Q/(k*gamma) + c1;
f4 = Q/(k*gamma*gamma);
den2 = (f1+f2)^2;
f3_h = (f1_h*(f1+f2) - f1*(f1_h+f2_h))/den2;
c1_h = f4*gamma_h*f3 - (Q/(k*gamma))*f3_h;
c2_h = -f4*gamma_h + c1_h;
c1_Q = -(1/(k*gamma))*f3;
c2_Q = (1/(k*gamma)) + c1_Q;

uvals = c1*exp(-gamma*xvals) + c2*exp(gamma*xvals) + u_amb;
uvals_data = c1*exp(-gamma*xdata) + c2*exp(gamma*xdata) + u_amb;
uvals_Q_data = c1_Q*exp(-gamma*xdata) + c2_Q*exp(gamma*xdata);
uvals_h_data = c1_h*exp(-gamma*xdata) + c2_h*exp(gamma*xdata) + gamma_h*xdata.*(-c1*exp(-gamma*xdata) + c2*exp(gamma*xdata));

res = data - uvals_data;
sens_mat = [uvals_Q_data; uvals_h_data];
sigma2 = (1/(n-p))*res*res';
V = sigma2*inv(sens_mat*sens_mat');



%%%3. Set MCMC parameters%%%
N = 1e+5;                        %number of iterations%
R = chol(V);
q_old = [Q;h];
SS_old = res*res';
n0 = 0.001;
sigma02 = sigma2;
aval = 0.5*(n0 + 15);
bval = 0.5*(n0*sigma02 + SS_old);
sigma2 = 1/gamrnd(aval,1/bval);
accept = 0;



%%%4. Run the Metropolis algorithm for N iterations.%%%
for i = 1:N
  z = randn(2,1); 
  q_new = q_old + R'*z;
  Q = q_new(1,1);
  h = q_new(2,1);
  gamma = sqrt(2*(a+b)*h/(a*b*k));
  f1 = exp(gamma*L)*(h + k*gamma);
  f2 = exp(-gamma*L)*(h - k*gamma);
  f3 = f1/(f2 + f1);
  c1 = -Q*f3/(k*gamma);
  c2 = Q/(k*gamma) + c1;
  uvals_data = c1*exp(-gamma*xdata) + c2*exp(gamma*xdata) + u_amb;
  res = data - uvals_data;
  SS_new = res*res';
  u_alpha = rand(1);
  term = exp(-.5*(SS_new-SS_old)/sigma2);
  alpha = min(1,term);
    if u_alpha < alpha
      Q_MCMC(:,i) = [Q; h];
      q_old = q_new;
      SS_old = SS_new;
      accept = accept + 1;
    else
      Q_MCMC(:,i) = q_old;
    end
  Sigma2(i) = sigma2;
  bval = 0.5*(n0*sigma02 + SS_old);
  sigma2 = 1/gamrnd(aval,1/bval);
end

Qvals = Q_MCMC(1,:);
hvals = Q_MCMC(2,:);

After running the code, we see in the figures below the chains and marginal densities for our parameters of interest.

In comparison, the chains and marginal densities that result from employing the Delayed Rejection Adaptive Metropolis (DRAM) Algorithm with 100,000 iterations are shown in the below figures.

We arrive at these plots with the following (this requires a function kde.m):

%%%5. Use kde to construct densities for Q and h.%%%
range_Q = max(Qvals) - min(Qvals);
range_h = max(hvals) - min(hvals);
Q_min = min(Qvals)-range_Q/10;
Q_max = max(Qvals)+range_Q/10;
h_min = min(hvals)-range_h/10;
h_max = max(hvals)+range_h/10;
[bandwidth_Q,density_Q,Qmesh,cdf_Q]=kde(Qvals);
[bandwidth_h,density_h,hmesh,cdf_h]=kde(hvals);
accept/N



%%%6. Plot solutions%%%
%Chain For \Phi Parameter%
figure(1)
plot(Qvals,'-','linewidth',2)
set(gca,'Fontsize',[22]);
axis([0 N -11 -9])           %L, R, D, U%
xlabel('Chain Iteration')
ylabel('Parameter \Phi')

%Chain For h parameter%
figure(2)
plot(hvals,'-','linewidth',2)
set(gca,'Fontsize',[22]);
axis([0 N 1.35e-3 1.5e-3])
xlabel('Chain Iteration')
ylabel('Parameter h')

%Error Variance%
figure(3)
plot(Sigma2)
set(gca,'Fontsize',[22]);
title('Measurement Error Variance \sigma^2')

%KDE for \Phi%
figure(4)
plot(Qmesh,density_Q,'k-','linewidth',3)
axis([-11 -9 0 5])
set(gca,'Fontsize',[22]);
xlabel('Parameter \Phi')
ylabel('PDF')

%KDE for h%
figure(5)
plot(hmesh,density_h,'k-','linewidth',3)
axis([1.3e-3 1.6e-3 0 4e4])
set(gca,'Fontsize',[22]);
xlabel('Parameter h')
ylabel('PDF')

figure(6)
scatter(Qvals,hvals)
box on
axis([-11 -9 1.3e-3 1.6e-3])
set(gca,'Fontsize',[22]);
xlabel('Parameter \Phi')
ylabel('Parameter h')

Both the Metropolis Algorithm and DRAM Algorithm show negative correlations between \(\Phi\) and \(h\), as can be seen from the sample points in the below figures.

The table below shows the estimates and standard deviations for our parameters of interest based on the three methods we have used so far (frequentist previously, the two Bayesian algorithms discussed above). Note that all three methods largely agree with one another.