%% Matlab Routine to Perform Bayesian Parameter Estimation
% Coded by Reuel Smith 2015-2017
% v. MATLAB R2015b through 2017a
% ========================================================================
% Example Problem 5.2
% Consider the following example involving a test of a particular metal
% alloy. The sliding test was conducted over a range of different applied
% weights in order to study the effect of weight and to gain a better
% understanding of the wear mechanism. The accelerated degradation test
% data is presented in Table 5.2 and plotted in log-log axes (Figure 5.8).
% Estimate the parameters of power degradation model using the MLE method,
% assuming no measurement errors. Repeat the problem assuming 10%
% multiplicative measurement error.
%
% ========================================================================
% Weight Cycles (hundreds)
% ==================================================================
% Unit 2 5 10 20 50 100 200 500
% ========================================================================
% 10g 1 3.2 4.1 4.5 4.7 5.8 6.8 7.7 9.6
% 2 2.7 3.4 3.8 3.9 5.4 5.7 6.3 8.4
% 3 2.1 2.7 3.1 3.3 4.0 4.6 5.7 6.6
% 4 2.6 3.5 4.0 4.0 5.2 6.1 6.7 8.5
% 50g 5 7.5 7.8 8.2 10.6 12.6 13.3 12.9 14.8
% 6 7.5 8.1 9.8 10.9 14.8 16.1 17.3 20.2
% 7 7.0 8.9 9.4 11.1 12.4 13.5 16.7 17.3
% 8 7.8 8.9 10.0 11.5 13.7 16.2 16.2 21.0
% 100g 9 12.5 15.4 17.2 20.5 24.1 27.0 29.4 37.9
% 10 11.0 13.9 16.1 18.6 22.2 27.8 31.0 36.6
% 11 13.0 15.1 18.6 20.2 23.9 29.7 31.5 39.6
% 12 11.7 13.7 16.7 17.5 22.3 25.3 32.0 38.2
% ========================================================================
% Solution:
% The example problem computes the mean and sigma for all the unit data.
% This example will only compute the parameter distribution of one of the
% units based on a Bayesian parameter analysis of
% log(D(cycles)) = log(sigma1) + sigma2*log(D(cycles))
% This is expected to follow a lognormal distribution.
% ========================================================================
clc
clear
% ========================================================================
% Initialize data and constants
% ========================================================================
% Definitions of different variables used in this code
% x(1): sigma 1 parameter
% x(2): sigma 2 parameter
% x(3): standard deviation parameter
% Cycle data (in hundreds of cycles)
N = [2 5 10 20 50 100 200 500];
Nfull = linspace(1,500,1000);
% Degradation measure width of scar (in microns). This degradation measure
% is used as the unit of damage at cycle N in this problem. This example
% only performs the Bayesian analysis on munit1, but the other munit data
% are available for experimenting with the code. C
munit1 = [3.2 4.1 4.5 4.7 5.8 6.8 7.7 9.6];
% munit2 = [2.7 3.4 3.8 3.9 5.4 5.7 6.3 8.4];
% munit3 = [2.1 2.7 3.1 3.3 4.0 4.6 5.7 6.6];
% munit4 = [2.6 3.5 4.0 4.0 5.2 6.1 6.7 8.5];
% munit5 = [7.5 7.8 8.2 10.6 12.6 13.3 12.9 14.8];
% munit6 = [7.5 8.1 9.8 10.9 14.8 16.1 17.3 20.2];
% munit7 = [7.0 8.9 9.4 11.1 12.4 13.5 16.7 17.3];
% munit8 = [7.8 8.9 10.0 11.5 13.7 16.2 16.2 21.0];
% munit9 = [12.5 15.4 17.2 20.5 24.1 27.0 29.4 37.9];
% munit10 = [11.0 13.9 16.1 18.6 22.2 27.8 31.0 36.6];
% munit11 = [13.0 15.1 18.6 20.2 23.9 29.7 31.5 39.6];
% munit12 = [11.7 13.7 16.7 17.5 22.3 25.3 32.0 38.2];
% m = [munit1;munit2;munit3;munit4;munit5;munit6;...
% munit7;munit8;munit9;munit10;munit11;munit12];
% number of parameters in likelihood equation: theta1, theta2, sigma (SD)
n=3;
% Initialize parameter guess
% This estimate is based on a polynomial fit of the log-linear equation
% which gives the theta1 and theta2 estimate. In a degradation equation
% where the relation is not linear, these values will have to be estimated
% in an alternate way such as an initial MLE estimate. Parameter SD is
% given as 1.
pinit = polyfit(log(N),log(munit1),1);
init_param=[pinit(2) pinit(1) 1];
% sigma for the proposed PDF, an nxn identity matrix
prop_sig=eye(n);
% lognormal PDF used in the likelihood. Change the expression
% x(1)+x(2)*log(N) when using a different damage equation
pdf1 = @(x,q) (1./((x(3)).*q.*sqrt(2*pi))).*exp(-0.5.*((log(q)-x(1)-x(2)*log(N))./x(3)).^2);
% number of randomly generated samples
nsamples = 20000;
%this value can remain as it is even if model parameters or model change
K = 1000;
% This parameter controls the size of the new markov chain which omits M-1
% out of M values. This will curb the effect of autocorrelation.
M = 10;
% define proposal distribution (proppdf) and random number generator
% (proprnd). The standard deviation for each parameter in the proprnd is
% to be chosen such that the updated guess covers a reasonable range of
% estimates. These should be small if the initial estimate is a close one.
proppdf = @(x,y) mvnpdf(x,y,prop_sig);
proprnd = @(y) [normrnd(y(1),0.015),normrnd(y(2),0.04),normrnd(y(3),0.0211)];
% define the LikelihoodxPriors
% The prior distributions for theta 1 and theta 2 are assumed to be
% non-informative, therefore they are treated as uniform distributions
% between -100 and 100. The prior distribution for SD is also a
% non-informative prior UNIF(0,2). These priors multiplied by the product
% of the lognormal PDF of each data point form the LikelihoodxPriors term
% that undergoes Bayesian updating.
pdf = @(x) unifpdf(x(1),-100,100)*unifpdf(x(2),-100,100)*unifpdf(x(3),0,2)*prod(pdf1(x,munit1));
% MH-sampling routine
[result,accept] = mhsample(init_param,nsamples,'pdf',pdf,'proppdf',proppdf,'proprnd',proprnd,'burnin',K,'thin',M);
% Alternatively the BURN-IN, and THIN value may be left out for a faster,
% but less precise MH-sampling procedure
% [result,accept] = mhsample(init_param,nsamples,'pdf',pdf,'proppdf',proppdf,'proprnd',proprnd);
% Posterior fit of theta 1 assuming a normal distribution
[mean1_mean,mean1_SD]=normfit(result(:,1));
% Posterior fit of theta 2 assuming a normal distribution
[mean2_mean,mean2_SD]=normfit(result(:,2));
% Posterior fit of sigma assuming a lognormal distribution
[mean3_meanSD]=lognfit(result(:,3));
% Best fit line based on the resulting means of theta 1 and theta 2
fbest = exp(polyval([mean2_mean mean1_mean],log(Nfull)));
% Plot of best fit line based on MH-Sample results
figure(1)
loglog(N,munit1,'b*',Nfull,fbest,'b-.')
xlabel('Cycles (hundreds)')
ylabel('Weight Loss (micrograms)')
grid on
set(gcf,'color','w');
figure(2)
plot(N,munit1,'b*',Nfull,fbest,'b-.')
xlabel('Cycles (hundreds)')
ylabel('Weight Loss (micrograms)')
grid on
set(gcf,'color','w');
figure(3)
subplot(1,3,1)
hist(result(:,1),100)
xlabel('ln(\theta_1)')
subplot(1,3,2)
hist(result(:,2),100)
xlabel('\theta_2')
subplot(1,3,3)
hist(result(:,3),100)
xlabel('\sigma')
set(gcf,'color','w');
Top