clear all

The following lines of code add a path to all the Bayesian statistical functions.

% Add MCMC codes to path
current_directory = cd;
addpath(strcat(current_directory,'/mcmcstat'));

Define a function that calculates the sum of squares of error between
the DFT energy and the approximation based on the Landau continuum
energy function.

%Define model for parameter estimation
% SS function and Parameters
ssfun = @coupled_func_Bayesian;

model.ssfun = ssfun;

Load the DFT calculations to fit continuum scale Landau energy.  A
scaling parameter is introduced (1e3) that converts the energy from GPa
units to MPa units.  This scaling is applied on the initial parameter
estimates and their bounds if the user wants to use a different set of
units in fitting Landau energy parameters.

load DFT_cubic
model.N  = length(Pnum); %length of the different polarization and energy states

scaling = 1e3; %all DFT energy density is in units of MPa using the scaling 1e3

%cubic volume
V0 = 4.1476034E+02; %volume of DFT lattice in Bohr^3
V1 = V0*(5.29177249e-11)^3; %conversion of volume to m^3
deltaE = (E_tot - E_tot(1))/V1*1.60217656535e-19/1e9; %Giga-Joules/m^3 conversion from eV/m^3

%input energy paramters from DFT
data.ydata_E_tot = deltaE*scaling;
data.xdata = Pnum; %polarization computed from DFT using the Berry phase approach
%define initial guess for Landau parameters
A = -0.805*scaling;
B = 3.32*scaling;
tmin = [A
        B];

%Landau energy model parameter range and initial guess
params = {
    {'A', tmin(1), -5e6*scaling, 5e6*scaling}
    {'B',  tmin(2), -5e3*scaling,5e8*scaling}
    };

% Example of output from this simulation
% MCMC statistics, nsimu = 2000
%
%                  mean         std      MC_err         tau      geweke
% ---------------------------------------------------------------------
%          A   -0.80461    0.015049  0.00080518      6.8223      0.9996
%          B     3.3196    0.037967   0.0016675      6.3636      0.9997
% ---------------------------------------------------------------------



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%check model initial output
%----------------------------------------------------------------
% THE NEXT LINES OF CODE CHECK THE PARAMETERS PRIOR TO RUNNING A PARAMETER
% ESTIMATION USING MCMC

[Sigma] = coupled_model_Bayesian(tmin,data.xdata);
ss = coupled_func_Bayesian(tmin,data);
SS = ss/length(data.xdata);  %here is the cost calculation

figure(1)
axes1 = axes('FontSize',20);
plot(Pnum,data.ydata_E_tot,'bo--','Linewidth',3)
hold on
plot(Pnum,Sigma,'r-','Linewidth',3)
hold off
ylabel('Energy (MJ/m^3)','FontSize',20)
xlabel('Polarization (C/m^2)','FontSize',20)
legend('Total energy (DFT)','Energy (continuum)')
legend1 = legend(axes1,'show');
set(legend1,'EdgeColor',[1 1 1],'YColor',[1 1 1],'XColor',[1 1 1]);

variance = sum((data.ydata_E_tot-Sigma').^2)/(length(deltaE-1));

%***************************************************************
% BAYESIAN ANALYSIS
model.S20 = variance;%estimate of variance based on prior simulation
model.N0  = 1;
model.N  = length(Pnum);
options.updatesigma=1;

options.nsimu = 500; %run an initial small set of iterations prior to larger run
[results, chain, s2chain] = mcmcrun(model,data,params,options);

options.nsimu = 2000; %run longer simulations using results from smaller run as
                      %initial conditions
[results, chain, s2chain] = mcmcrun(model,data,params,options,results);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot results
figure(2)
mcmcplot(chain,[],results,'denspanel',2);
figure(3)
mcmcplot(chain,[],results,'pairs');

figure(4); clf
mcmcplot(chain,[],results.names,'chainpanel')
xlabel('Iterations','Fontsize',14)
ylabel('Parameter value','Fontsize',14)

chainstats(chain,results) %summary of parameter means, standard deviations, etc.

%Define a model function that calls the chains for error propagation
modelfun1 = @(d,th)coupled_model_Bayesian(th,d); % NOTE: the order in which d and th appear is important

%The following lines of code propagate the posterior densities of
%parameters to determine prediction and credible intervals of the Landau
%energy with respect to the DFT energy in light of the Bayesian statistics
nsample = 500;
out = mcmcpred(results,chain,s2chain,data.xdata,modelfun1,nsample);

figure(5)
modelout = mcmcpredplot(out); %,data
hold on
plot(data.xdata,data.ydata_E_tot,'bo-','linewidth',3)
hold off
xlabel('P_{3} (C/m^2)','Fontsize',14);
ylabel('Energy (MPa)','FontSize',20)
legend('95/% Prediction Interval','95% Confidence Interval','Continuum','DFT','Location','Best')

Setting nbatch to 1
Sampling these parameters:
name   start [min,max] N(mu,s^2)
A: -805 [-5e+09,5e+09] N(0,Inf)
B: 3320 [-5e+06,5e+11] N(0,Inf)
Setting nbatch to 1
Using values from the previous run
Sampling these parameters:
name   start [min,max] N(mu,s^2)
A: -803.278 [-5e+09,5e+09] N(0,Inf)
B: 3298.9 [-5e+06,5e+11] N(0,Inf)
MCMC statistics, nsimu = 2000

                 mean         std      MC_err         tau      geweke
---------------------------------------------------------------------
         A     -803.9      17.255       0.969      7.9091     0.99935
         B     3317.9      42.415      2.3866      7.6127     0.99825
---------------------------------------------------------------------

Published with MATLAB® R2014a