Bayesian updating for a simple concrete creep model

Hello,

Here is my first try with bayesian updating in openturns. This is open for discussion !

The model is a concrete creep model (called beton_burger_fp in Code_Aster) for an uniaxial creep test during 60 years. Data are available for 10 years.

6 inputs are k_rs, k_rd, eta_rs, eta_rd, eta_id, kappa.

Outputs are:

  1. a priori versus a posteriori inputs distribution
  2. a priori and a posteriori mean response and 95%-quantile of the response.

The main script is actu_bayes.py. fluage.data and StepB.py are used by actu_bayes.py.

fluage.data:

 0.              0.
2.46E-09       -1.71E-19
3.74E+03       2.20E-05
1.95E+04       1.15E-04
1.29E+05       3.54E-05
1.34E+05       6.47E-05
2.69E+05       1.40E-04
1.01E+06       1.99E-04
1.88E+06       2.35E-04
2.12E+06       2.43E-04
2.37E+06       2.54E-04
2.74E+06       2.60E-04
3.11E+06       2.78E-04
4.46E+06       2.99E-04
5.32E+06       3.07E-04
5.94E+06       3.18E-04
6.18E+06       3.37E-04
7.53E+06       3.44E-04
8.15E+06       3.50E-04
8.64E+06       3.56E-04
9.62E+06       3.65E-04
1.07E+07       3.83E-04
1.21E+07       4.03E-04
1.32E+07       4.26E-04
1.53E+07       4.23E-04
1.60E+07       4.38E-04
1.85E+07       4.40E-04
1.99E+07       4.50E-04
2.26E+07       4.78E-04
2.55E+07       4.85E-04
2.59E+07       5.07E-04
2.94E+07       5.07E-04
3.34E+07       5.27E-04
3.44E+07       5.39E-04
3.72E+07       5.57E-04
4.07E+07       5.71E-04
4.90E+07       6.08E-04
7.55E+07       6.92E-04
9.46E+07       7.60E-04
1.14E+08       8.25E-04
1.32E+08       8.58E-04
1.61E+08       9.87E-04
1.99E+08       1.12E-03
2.27E+08       1.21E-03
2.55E+08       1.31E-03
2.84E+08       1.41E-03
3.12E+08       1.51E-03

StepB.py:

  # -*- coding: utf-8 -*-

from openturns import *

import numpy as np

dim_in  = 7        # Nombre de variables d'entree
dim_out = nb_data+nb_disc  # Nombre de variables de sortie

# Fonction de transfert

class myfunction(OpenTURNSPythonFunction):
        def __init__(self):
                OpenTURNSPythonFunction.__init__(self, dim_in, dim_out)
        def _exec(self, X):     
                        
                v_KRS   = X[0];
                v_KRD   = X[1];
                v_ETARS = X[2];
                v_ETARD = X[3];
                v_ETAID = X[4];
                v_KAPPA = X[5];

                nu = 0.2; #coefficient de poisson elastique

                T_FIN = 1.892E9; #60 ans

                TEMPS = np.zeros(dim_out)
  
                for i in range(0,len(temps)):
                  TEMPS[i] = np.array(temps[i])

                T_FIN_data = TEMPS[nb_data-1]

                SIGMA_0 = 12E6;

                ETA_IS = (1+nu)/(1-2*nu)*v_ETAID
                B = v_KAPPA/(np.sqrt(1.+2.*nu**2));
                C = SIGMA_0*np.sqrt(1.+2.*nu**2)/(v_KAPPA*(1.+nu)*v_ETAID);
                C2 = SIGMA_0*np.sqrt(1.+2.*nu**2)/(v_KAPPA*(1.-2.*nu)*ETA_IS);
                D = -v_KRD/v_ETARD;
                D2 = -v_KRS/v_ETARS;

                EPSA = [];
                EPS1 = [];

                for i in range(1,nb_disc+1):
                        TEMPS[nb_data+i-1] = T_FIN_data+i*(T_FIN-T_FIN_data)/nb_disc
                for i in range(dim_out):
                        EPSA.append(-SIGMA_0/3.*(1./v_KRS*(1.-np.exp(D2*TEMPS[i]))+2./v_KRD*(1.-np.exp(D*TEMPS[i])))-B*np.log(1.+(3.*SIGMA_0*np.sqrt(1.+2.*nu**2)*TEMPS[i])/(v_KAPPA*(ETA_IS*(1.-2.*nu)+2.*v_ETAID*(1.+nu)))));
                        EPS1.append(-EPSA[i]);


                return EPS1

model = NumericalMathFunction(myfunction())

actu_bayes.py

# -*- coding: utf-8 -*-
from openturns import *
import numpy as np
import pylab as pl
from scipy.stats import norm
import time
import scipy.optimize as so
import scipy.stats as st

###################
# Parametres civaux
###################

#Proprietes elastiques
nu = 0.2; #coefficient de poisson elastique
sigma_0 = 12E6;

#########
# Instants des donnees experimentales
#########

data       = open('fluage.data','r')

liste = []
liste = data.readlines()
nb_data = len(liste)
data.close()

temps = NumericalSample(nb_data,1)

i=0
for j in range(len(liste)):
  ligne1 = liste[j].split()
  temps[i,0]=float(ligne1[0])
  i = i+1

######
# Definition du modele
######

def buildModel(temps):

  # Express eta_0_is in terms of the other parameters
  eta_0_is = "eta_0_id * (" + str((nu + 1.0) / (1.0 - 2.0 * nu)) + ")"
  # Convert all the fixed numerical values into strings
  # Protect them with parenthesis to avoid troubles with negative values
  sigma_0_str = str(sigma_0)
  t_str = "(" + str(temps) + ")"
  factor = np.sqrt(1.0 + 2.0 * nu**2)
  term0 = "(" + str(factor) + ")"
  term1 = "(" + str(3.0 * sigma_0 * factor) + ")"
  term2 = "(" + str(1.0 - 2.0 * nu) + ")"
  term3 = "(" + str(1.0 + nu) + ")"
  # The expression, at last...
  return NumericalMathFunction(["k_rs", "k_rd", "eta_rs", "eta_rd", "eta_0_id", "kappa", "log_sigma"], [sigma_0_str + " / 3 * (1.0 / k_rs * (1.0 - exp(-k_rs / eta_rs * " + t_str + ")) + 2 / k_rd * (1.0 - exp(-k_rd / eta_rd * " + t_str + "))) + kappa / " + term0 + " * log(1.0 + " + term1 + " / (kappa * (" + eta_0_is + " * " + term2 + " + 2.0 * eta_0_id * " + term3 + ")) * " + t_str + ")", "exp(log_sigma)"])

# Version Numpy du modele 

T = np.array(temps).ravel()

def modeleDeformation(theta):
    k_rs,k_rd,eta_rs,eta_rd,eta_id,kappa, log_sigma = theta

    eta_is = (1+nu)/(1-2*nu)*eta_id
    B   = kappa/(np.sqrt(1+2*nu*nu))
    C   = sigma_0*np.sqrt(1+2*nu*nu)/(kappa*(1+nu)*eta_id)
    C2  = sigma_0*np.sqrt(1+2*nu*nu)/(kappa*(1-2*nu)*eta_is)
    D   = -k_rd/eta_rd
    D2  = -k_rs/eta_rs

    return np.hstack((sigma_0/3*(1/k_rs*(1-np.exp(D2*T))+2/k_rd*(1-np.exp(D*T))) \
            + B*np.log(1+(3*sigma_0*np.sqrt(1+2*nu**2)*T)/(kappa*(eta_is*(1-2*nu)+2*eta_id*(1+nu)))),\
            np.exp(log_sigma)))

###
# Mesures de deformation
###

yobs= [ 0., -1.71E-19, 2.20E-05, 1.15E-04, 3.54E-05, 6.47E-05, 1.40E-04, 1.99E-04, \
 2.35E-04, 2.43E-04, 2.54E-04, 2.60E-04, 2.78E-04, 2.99E-04, 3.07E-04, 3.18E-04, \
 3.37E-04, 3.44E-04, 3.50E-04, 3.56E-04, 3.65E-04, 3.83E-04, 4.03E-04, 4.26E-04, \
 4.23E-04, 4.38E-04, 4.40E-04, 4.50E-04, 4.78E-04, 4.85E-04, 5.07E-04, 5.07E-04, \
 5.27E-04, 5.39E-04, 5.57E-04, 5.71E-04, 6.08E-04, 6.92E-04, 7.60E-04, 8.25E-04, \
 8.58E-04, 9.87E-04, 1.12E-03, 1.21E-03, 1.31E-03, 1.41E-03, 1.51E-03]

obsSize = len(yobs)

yobs = np.atleast_2d(yobs).T
yobs = NumericalSample(yobs)

modelCollection=NumericalMathFunctionCollection()

for iaux in range(nb_data):
  modelCollection.add(buildModel(temps[iaux,0]))

model = NumericalMathFunction(modelCollection)

conditional = Normal()

# reference
Moyenne = [2.26E10, 3.46E10, 1.97E17, 1.E17, 2.51E18, 2.51E-3]

# on donne plus ou moins de poids au modele : 1/sigma^2
# variable tres importante
StdDev = NumericalPoint([elt for elt in NumericalPoint(np.hstack(([0.15*elt for elt in Moyenne],-np.log(.00025))))])
CorrMatrix = CorrelationMatrix(7)

for iaux in range(7):
  CorrMatrix[iaux,iaux] = StdDev[iaux]*StdDev[iaux]

Prior = Normal(NumericalPoint(np.hstack((Moyenne,.00025))),CorrMatrix)

Theta0 = Prior.getMean()

print "Moyenne prior " , NumericalPoint(np.hstack((Moyenne,.00025)))
print "Prior.getMean : " , Prior.getMean()


proposal = DistributionCollection()
for iaux in range(len(StdDev)):
    width = 0.15*StdDev[iaux]
    proposal.add(Uniform(-width,width))

RWMHsampler = RandomWalkMetropolisHastings(Prior,conditional,model,yobs,Theta0,proposal)

strategy = CalibrationStrategyCollection()
aStrategy = CalibrationStrategy(Interval(0.,1.))
for iaux in range(7):
  strategy.add(aStrategy)

RWMHsampler.setCalibrationStrategyPerComponent(strategy)

RWMHsampler.setVerbose(True)
RWMHsampler.setThinning(1)
RWMHsampler.setBurnIn(100)

sampleSize = 200

t1 = time.time()
sample = RWMHsampler.getSample(sampleSize)
t2 = time.time()
print 'acceptance rate = ',
print RWMHsampler.getAcceptanceRate()


myRandomVector = PosteriorRandomVector(RWMHsampler)
# AUGMENTER A 100000 pour avoir des distributions tres plates
n = 5000

posterior_sample = myRandomVector.getSample(n)
np_posterior_sample = np.array(posterior_sample)
#np_posterior_sample = np.array(sample)

post_mean = np_posterior_sample.mean(axis=0)

#Posterior plot
nsample = 1000

pl.clf()
pl.figure(999)
for i in xrange(nsample):
    bb = pl.plot(T/3600./24./365.,modeleDeformation(np_posterior_sample[np.random.randint(n)])[:-1], 'r:')

cc = pl.plot(T/3600./24./365.,modeleDeformation(post_mean)[:-1], 'k', lw=2)
aa = pl.scatter(T/3600./24./365.,yobs, s=20, c='g',zorder=10)

pl.legend((aa, bb, cc), ("obs", "Loi a posteriori", "Moyenne a posteriori"), loc="upper left")
pl.xlabel("Temps (annees)")
pl.ylabel("Deformation differee")
pl.savefig("AjustementMoyenneAPosteriori.pdf")
pl.xlim(0,3.)
pl.ylim(0,0.001)
pl.savefig("AjustementMoyenneAPosterioriZoom.pdf")


pl.figure()
for i in range(7):
    #pl.figure(50+i)
    pl.subplot(7,1,i+1)
    pl.plot(range(n), np_posterior_sample[:,i])

filename = "cv_posterior_var.pdf"
pl.savefig(filename)

pl.clf()

t3 = time.time()

kernel = KernelSmoothing()

# Prior input PDFs (for comparison)
for i in range(7):
  priorPDF = Prior.getMarginal(i).drawPDF()
  priorPDFdrawable = priorPDF.getDrawable(0)
  posteriordrawable = kernel.build(posterior_sample.getMarginal(i)).drawPDF().getDrawable(0) 
  posteriordrawable.setLegendName('posterior')
  posteriordrawable.setColor('blue')
  priorPDF.add(posteriordrawable)
  repe = "prior_vs_posterior_var_entree_"+str(i+1)+"_sampleSize"+str(sampleSize)+"_n_"+str(n)
  if not(repe in os.listdir('.')) :
    os.mkdir(repe)
  os.chdir(repe)
  filename = "prior_vs_posterior_var_entree_"+str(i+1)+"_sampleSize"+str(sampleSize)+"_n_"+str(n)
  priorPDF.draw(filename)
  os.chdir('../.')

np.savetxt('echantillon_'+str(n),np_posterior_sample)
print 'calcul de la loi a posteriori : ',t2-t1
print 'calcul d un echantillon a posteriori : ',t3-t2

print "LogVraisemblance(Theta0) = ", RWMHsampler.computeLogLikelihood(Theta0)
print "LogVraisemblance(10*Theta0) = ", RWMHsampler.computeLogLikelihood(10*Theta0)

pl.clf()
pl.figure(7)
####
# Impression des donnees
####
data       = open('fluage.data','r')
temps_data = []
eps1_data  = []

liste = []
liste = data.readlines()
data.close()

for j in range(1,len(liste)):
  ligne1 = liste[j].split()
  temps_data.append(float(ligne1[0]))
  eps1_data.append(float(ligne1[1]))
for ind in range(len(temps_data)):
  temps_data[ind] = temps_data[ind]/(3600.*24.*365.)

pl.plot(temps_data, eps1_data,'bo' ,label="Donnees essais exp.")
pl.legend(bbox_to_anchor=(.05, 1), loc=2, borderaxespad=0.)
pl.xlabel("Temps (annees)")
pl.ylabel("Deformation differee")
pl.title("Comparaison actu bayes et essais")

###
# Impression du posterior et de l intervalle de confiance a 95 %
####
# discretisation temporelle au dela des instants pour lesquels on dispose de
# donnees experimentales
nb_disc = 600

execfile('StepB.py')

sizeX = n

SampleMC     = model(posterior_sample)
priorSample  = Prior.getSample(sizeX)
sample_model_prior = model(priorSample)

meanMC       = SampleMC.computeMean()  
covarianceMC = SampleMC.computeCovariance()
meanMCprior       = sample_model_prior.computeMean()  
covarianceMCprior = sample_model_prior.computeCovariance()

# calcul de l'intervalle de confiance a 95% de l'estimateur meanMC
borneinf_meanMC=[]
bornesup_meanMC=[]

alpha = 0.95

lower_bnds=np.array(SampleMC.computeQuantilePerComponent((1-alpha)/2.) )
upper_bnds=np.array(SampleMC.computeQuantilePerComponent((1+alpha)/2.) )
lower_bnds_pr=np.array(sample_model_prior.computeQuantilePerComponent((1-alpha)/2.) )
upper_bnds_pr=np.array(sample_model_prior.computeQuantilePerComponent((1+alpha)/2.) )

for i in range(dim_out):
  borneinf_meanMC.append(meanMC[i]-1.96*np.sqrt(covarianceMC[i, i])/np.sqrt(sizeX))
  bornesup_meanMC.append(meanMC[i]+1.96*np.sqrt(covarianceMC[i, i])/np.sqrt(sizeX))

########
# On ramene en annees
########

TEMPS   = [];

# Instant final 60 ans
T_FIN = 1.892E9;
TEMPS = np.zeros(dim_out)
for i in range(0,len(temps)):
    TEMPS[i] = np.array(temps[i])
T_FIN_data = TEMPS[nb_data-1]
for i in range(1,nb_disc+1):
    TEMPS[nb_data+i-1] = T_FIN_data+i*(T_FIN-T_FIN_data)/nb_disc

for ind in range(len(TEMPS)) :
  TEMPS[ind] = TEMPS[ind]/(3600.*24.*365.)

pl.clf()
pl.figure(100)
pl.xlabel("Temps (annees)")
pl.ylabel("Deformation differee")
pl.title("Comparaison actu bayes et essais")
pl.plot(temps_data, eps1_data,'bo' ,label="Donnees essais exp.")
pl.plot(TEMPS, meanMC,'g', label="Moyenne MC posterior")
pl.plot(TEMPS, meanMCprior,'k', label="Moyenne MC prior")
pl.plot(TEMPS,lower_bnds,'b',label="quantile 2.5 posterior")
pl.plot(TEMPS,upper_bnds,'b',label="quantile 97.5 posterior")
pl.plot(TEMPS,lower_bnds_pr,'r',label="quantile 2.5 prior")
pl.plot(TEMPS,upper_bnds_pr,'r',label="quantile 97.5 prior")
filename = "estimation_cv_MC_taille_"+str(sizeX)
pl.legend(loc='best')
pl.savefig(filename)
pl.xlim(0,3.)
pl.ylim(0,0.001)
filename="estimation_cv_MC_taille_restreint_"+str(sizeX)
pl.savefig(filename)

for ind in range(len(TEMPS)) :
  TEMPS[ind] = TEMPS[ind]*(3600.*24.*365.)

print 'moyenne de l echantillon posteriori = ',posterior_sample.computeMean()
print 'ecart type de l echantillon posteriori = ',posterior_sample.computeStandardDeviationPerComponent()

for i in range(6):
  print 'coefficient de variation variable '+str(i+1)+' : ',posterior_sample.computeStandardDeviationPerComponent()[i]/posterior_sample.computeMean()[i]

Comments

No comments.