The supplemental materials consist of two MATLAB scripts, MATLAB_Code_S1.m and MATLAB_Code_S2.m. They are provided below, as well as in separate .m MATLAB files. The two files are subroutines that provide calculations used in the development of the site-specific seismological model presented in the main article.
MATLAB_Code_S1.m calculates the basis functions for the spectral density function gj(f) as shown in equation (5) in the main article.
MATLAB_Code_S2.m performs the update of the unknown parameters (Θ1, Θ2, and Θ3) as part of the stochastic surrogate model for the spectral density function shown in equation (6). The update is done according to the development described in the main article in the Bayesian Analysis section and following equations (11)–(15).
Download: MATLAB_Code_S1.m.zip [Zipped MATLAB Source Code; 1 KB]. This is a downloadable version of the code shown below.
% =====================================================================
% Surrogate Model for the SBM Spectral Densities
% -----------------------------------------------
% by Alin Radu, 05/13/2014
% e-mail: acr99@cornell.edu
% Cornell University
% Ithaca, NY 14850
% =====================================================================
% Description: It provides the orthogonal decomposition of the SBM spectral
% density. The modes of the decomposition are obtained through singular
% value decomposition and they are grouped in three deterministic functions
% phi1, phi2, phi3. The spectral density function for a fixed (m,r) can be
% written as:
% gw(w;m,r)=phi1(w;m,r)+phi2(w;m,r)+phi3(w;m,r)
% =====================================================================
% INPUT
% The input is stored in the SBM_psd.mat file which contains:
% cells = matrix of dimentions (N × 2), where N is the total number of
% cells used in the seismic activity matrix, i.e., the total number
% of pairs (m,r). The first column of cells is the moment magnitude
% (m) and the second column is the source-to-site distances (r).
% w = vector of discrete frequency points
% gw = the spectral densities for each pair of (m,r) as obtained from the
% seismological model based on the specific barrier model put forth in
% a seismological model by Halldorson & Papageorgiu (2005). Size of
% gw is (length(w) × N).
% ------------------------------
% OUTPUT
% phi1, phi2, phi3 = three modes which compose the spectral density, i.e.,
% gw(w;m,r)=phi1(w;m,r)+phi2(w;m,r)+phi3(w;m,r)
% =====================================================================
function [phi1,phi2,phi3]=MATLAB_Code_S1(w,gw,cells)
[uk,sk,vk]=svd(gw);
phi1=zeros(size(sk));
phi2=zeros(size(sk));
phi3=zeros(size(sk));
for icell=1:length(cells)
phi1(:,icell)=uk(:,1)*vk(icell,1)*sk(1,1);
phi2(:,icell)=uk(:,2)*vk(icell,2)*sk(2,2);
for i=3:length(uk)
phi3(:,icell)=phi3(:,icell)+sk(i,i)*uk(:,i)*vk(icell,i);
end
end
end
Download: MATLAB_Code_S2.m.zip [Zipped MATLAB Source Code; 1 KB]. This is a downloadable version of the code shown below.
% =====================================================================
% Statistical Update of the SBM Spectral Densities
% -----------------------------------------------
% by Alin Radu, 05/13/2014
% e-mail: acr99@cornell.edu
% Cornell University
% Ithaca, NY 14850
% =====================================================================
% Description: It provides the updated version of the porbability density
% function for the random parameters (theta1, theta2, theta3) of the
% stochastic surrogate model for the spectral density:
% ~gw(w;m,r)=theta1*phi1(w;m,r)+theta2*phi2(w;m,r)+theta3*phi3(w;m,r)
% =====================================================================
% INPUT
% w = vector of discrete frequency points
% phi1, phi2, phi3 = modes for the SBM power spectral density g(w;m,r) in
% the desired cell of the seismic activity matrix for
% the characteristics (m,r)_j. Sizes of phi1, phi2, phi3
% are the same as the size of w.
% theta1, theta2, theta3 = vectors with the ranges for the unknown
% parameters.
% prior = joint prior probability density for the unknown parameters Theta,
% where Theta=(theta1,theta2,theta3).
% Size of prior is (length(theta1),length(theta2),length(theta3)).
% t = time vector for the ground motion sample used to update Theta.
% x = ground motion sample for (m,r)_j. Sizes of t and x are identical.
% ------------------------------
% OUTPUT
% posterior = joint posterior probability density for the unknown
% parameters Theta after the update to the ground motion record
% x(t).
% =====================================================================
function [posterior]=MATLAB_Code_S2(w,phi1,phi2,phi3,prior,t,x,theta1,theta2,theta3)
nth1=length(theta1);
nth2=length(theta2);
nth3=length(theta3);
l1=zeros(size(prior));
l2=zeros(size(prior));
dt=pi/w(end);
nt=ceil(t(end)/dt);
% tfin=nt*dt;
n=nt+1;
c=zeros(n);
count=0;
for i1=1:nth1
for i2=1:nth2
for i3=1:nth3
count=count+1;
gg=phi1*theta1(i1)+phi2*theta2(i2)+phi3*theta3(i3);
for i=1:n
for j=i:n
dt=t(i)-t(j);
c(i,j)=trapz(w,gg’.*cos(w*dt));
c(j,i)=c(i,j);
end
end
[vect,vals]=eig(c);
vals=diag(vals);
l1(i1,i2,i3)=−0.5*(x*(c\eye(n))*x’);
l2(i1,i2,i3)=−0.5*sum(log(vals(vals>0)));
end
end
end
like=l1+l2;
like=like.*prior;
posterior=like/sum(sum(sum(like)));
end
[ Back ]