The following is the code for my model. To use these code, need to install ANSYS as a outer calculator

```
import numpy as np
import pytensor.tensor as pt
calls = 0
def fem(mapdl, parameters):
mapdl.clear()
mapdl.prep7()
x = [1] * 20
x[0] = parameters[0] # 一层柱
x[5] = parameters[0] # 二层柱
x[10] = parameters[2] # 三层柱
x[15] = parameters[2] # 四层柱
y = [1] * 20
y[3] = parameters[1] # 一层梁
y[8] = parameters[1] # 二层梁
y[13] = parameters[3] # 三层梁
y[18] = parameters[3] # 四层梁
for i in range(20):
mapdl.et(i + 1, 'beam189')
mapdl.mp("EX", i + 1, x[i]*2.1e11) # N/mm2
mapdl.mp('DENS', i + 1, y[i]*7.8e3) # kg/mm3
mapdl.mp("PRXY", i + 1, 0.3) # Poisson's ratio
mapdl.sectype(i + 1, "BEAM", "RECT")
mapdl.secdata(0.065, 4e-3)
for i in range(5):
mapdl.k(3 * i + 1, -0.350, 0.350 * i, 0)
mapdl.k(3 * i + 2, 0, 0.350 * i, 0)
mapdl.k(3 * i + 3, 0.350, 0.350 * i, 0)
mapdl.k(16, 0.700, 0.700, 0)
mapdl.k(17, 0, 0.7700, 0)
for i in range(4):
mapdl.l(3 * i + 1, 3 * i + 4)
mapdl.l(3 * i + 2, 3 * i + 5)
mapdl.l(3 * i + 3, 3 * i + 6)
mapdl.l(3 * i + 4, 3 * i + 5)
mapdl.l(3 * i + 5, 3 * i + 6)
# 将材料、实常数、单元类型、、方向点、截面类型赋予此前选中的线模型
for i in range(4):
# 选中各层柱线,1一层柱，6二层柱，11三层柱，16四层柱
mapdl.lsel("S", "LINE", "", [5 * i + 1, 5 * i + 2, 5 * i + 3])
# 将材料、实常数、单元类型、、方向点、截面类型赋予此前选中的线模型
mapdl.latt(5 * i + 1, 5 * i + 1, 5 * i + 1, "", 16, 5 * i + 1)
# 选中各层梁线
mapdl.lsel("S", "LINE", "", [5 * i + 4, 5 * i + 5])
mapdl.latt(5 * i + 4, 5 * i + 4, 5 * i + 4, "", 17, 5 * i + 4)
# 网格划分
mapdl.allsel()
mapdl.lesize("all", "", "", 1)
mapdl.lmesh("all")
# 执行分析
mapdl.slashsolu()
mapdl.dk(1, "all")
mapdl.dk(2, "all")
mapdl.dk(3, "all")
mapdl.modal_analysis(nmode=6, freqb=1)
# ==============获取结果
disps = []
for i in range(3):
displacement = mapdl.result.nodal_displacement(i)
disp = displacement[1][:, 0]
disp = disp / np.linalg.norm(disp, ord=2)
disps.append(disp)
mapdl.post1()
freq = np.zeros(6)
for i in range(6):
mapdl.set(1, i+1)
freq[i] = mapdl.post_processing.freq
mapdl.finish()
return freq
```

this is the observed data I artificially created.

```
# make observed data, input=numpy, output=numpy list NS*6
def real(mapdl, NS, parameters):
# 定义各独立随机变量的参数
samples = np.zeros((NS, len(parameters))) # 创建一个数组来存储样本
# 对参数采样并计算均值
print("参数样本均值:")
for i in range(len(parameters)):
samples[:, i] = np.random.normal(loc=parameters[i]["mean"], scale=parameters[i]["std"], size=NS)
print(np.mean(samples[:, i]))
# 定义储存实测数据的list，并储存
freq_list = []
#disp1_list = []
#disp2_list = []
#disp3_list = []
for _NS in range(NS):
frequent = fem(mapdl, samples[_NS])
freq_list.append(frequent)
#disp1_list.append(modal1)
#disp2_list.append(modal2)
#disp3_list.append(modal3)
freq_array = np.array(freq_list)
return freq_array
```

```
# define a loglike function to process the data and the simulate
def my_loglike(mapdl, r, data, sigma):
simu_freq = fem_wrapper(mapdl, r)
loglike = 0
for _log in range(len(data)):
loglike += -(0.5 / sigma**2) * np.sum((data[_log] - simu_freq) ** 2)
return loglike
# define a pytensor Op for likelihood function
class LogLike(pt.Op):
itypes = [pt.dvector] # expects a vector of parameter values when called
otypes = [pt.dscalar] # outputs a single scalar value (the log likelihood)
def __init__(self, mapdl, loglike, data, sigma):
# add inputs as class attributes
self.mapdl = mapdl
self.likelihood = loglike # 传递似然函数
self.data = data # 观测到的数据
self.sigma = sigma # 似然函数的标准差
def perform(self, node, inputs, outputs):
# the method that is used when calling the Op
(parameters,) = inputs # this will contain my variables
# call the log-likelihood function
logl = self.likelihood(self.mapdl, parameters, self.data, self.sigma)
outputs[0][0] = np.array(logl) # output the log-likelihood
def fem_wrapper(mapdl, parameters):
global calls
calls += 1
return fem(mapdl, parameters)
def count():
global calls
print("The number of calls to FEM is", calls)```
```