Hi there! im trying to eastimate parameters of a structure , this is my first time to write py code.
i got a response of the structure like this
I knew how F(x）work and i want to estimate 10 parameters of the structure.
in this picture
Y=F(Sigma(i),noise)+ white gaussian noise
so I create a vitural model to estimate first.
I create Y=F(sigma=1,noise)+ white noise by matlab to estimate whether Sigma(i)=1, this is my code
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import pandas as pd
d_1 = pd.read_csv('d1.csv')#data of Y I create
noise = pd.read_csv('noise.csv')#force input
#get the data of structure
######--------model##################
with pm.model() as model:
#--------put in parameters----------------#
sigma1 = pm.Normal('sigma1', mu=0.9, sd=1)
sigma2 = pm.Normal('sigma2', mu=0.9, sd=1)
sigma3 = pm.Normal('sigma3', mu=0.9, sd=1)
sigma4 = pm.Normal('sigma4', mu=0.9, sd=1)
sigma5 = pm.Normal('sigma5', mu=0.9, sd=1)
sigma6 = pm.Normal('sigma6', mu=0.9, sd=1)
sigma7 = pm.Normal('sigma7', mu=0.9, sd=1)
sigma8 = pm.Normal('sigma8', mu=0.9, sd=1)
sigma9 = pm.Normal('sigma9', mu=0.9, sd=1)
sigma10 = pm.Normal('sigma10', mu=0.9, sd=1)
#---------------the F[x] to get mus----------------#
K = np.matrix([[sigma1+sigma2, -sigma2, 0, 0, 0, 0, 0, 0, 0, 0],
[-sigma1, sigma2+sigma3, -sigma3, 0, 0, 0, 0, 0, 0, 0],
[0, -sigma2, sigma3+sigma4, -sigma4, 0, 0, 0, 0, 0, 0],
[0, 0, -sigma3, sigma4+sigma5, -sigma5, 0, 0, 0, 0, 0],
[0, 0, 0, -sigma4, sigma5+sigma6, -sigma6, 0, 0, 0, 0],
[0, 0, 0, 0, -sigma5, sigma6+sigma7, -sigma7, 0, 0, 0],
[0, 0, 0, 0, 0, -sigma6, sigma7+sigma8, -sigma8, 0, 0],
[0, 0, 0, 0, 0, 0, -sigma7, sigma8+sigma9, -sigma9, 0],
[0, 0, 0, 0, 0, 0, 0, -sigma8, sigma9+sigma10, -sigma10],
[0, 0, 0, 0, 0, 0, 0, 0, -sigma9, sigma10]])
M = np.identity(10)
#mass matrix
mass = np.ones(10, 10, 1)
C = 0.05*M+0.02*K
#rayleigh
dt = 0.01
N = 1000
t = np.linspace(0.0, dt*(N-1), N)
#full time
#newmark method
alpha = 0.25
beta = 0.5
a0 = 1/alpha/dt/dt
a1 = beta/alpha/dt
a2 = 1/alpha/dt
a3 = 1/2/alpha-1
a4 = beta/alpha-1
a5 = dt/2*(beta/alpha-2)
a6 = dt*(1-beta)
a7 = dt*beta
g = 9.8
#get F(x) at t=0
deal = np.zeros(10, 1000)
vela = np.zeros(10, 1000)
areal = np.zeros(10, 1000)
#------- round to get real data -----------
for i in range(1, 999):
f = M*noise[i-1]*g
Ke = K+a0*M+a1*C
Fe = f+M*(a0*deal[:, i-1]+a2*vela[:, i-1]+a3*areal[:, i-1])+C*(a1*deal[:, i-1]+a4*vela[:, i-1]+a5*areal[:, i-1])
deal[:, i] = np.linalg.inv(Ke)*Fe
areal[:, i] = a0*(deal[:, i]-deal[:, i-1])-a2*vela[:, i-1]-a3*areal[:, i-1]
vela[:, i] = vela[:, i-1]+a6*areal[:, i-1]+a7*areal[:, i]
d1 = deal[0, :]
####### get mus already #################
mus = np.ones(10) * 0.9
cov = np.identity(10)
prior = pm.MvNormal('prior', mu=mus, cov=cov, shape=10)
y = pm.Normal('y', mus=d1, sd=1, observed=d_1)
with model:
trace = pm.sample(1000)
pm.traceplot(trace, ['sigma1', 'sigma2', 'sigma3', 'sigma4', 'sigma5', 'sigma6', 'sigma7', 'sigma8', 'sigma9'])
pm.traceplot(trace, ['sigma10'])
plt.show()
is this code right?
but i got warning in pycharm, how solve this ?
WARNING (theano.configdefaults): g++ not available, if using conda: `conda install m2w64-toolchain`
C:\Users\Icesp\PycharmProjects\chrispy\venv\lib\site-packages\theano\configdefaults.py:560: UserWarning: DeprecationWarning: there is no c++ compiler.This is deprecated and with Theano 0.11 a c++ compiler will be mandatory
warnings.warn("DeprecationWarning: there is no c++ compiler."
WARNING (theano.configdefaults): g++ not detected ! Theano will be unable to execute optimized C-implementations (for both CPU and GPU) and will default to Python implementations. Performance will be severely degraded. To remove this warning, set Theano flags cxx to an empty string.
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.