Hello, I am just in the initial stage of discovering pyMC. So, please forgive me if my question is really basic or not understandable. I have not found the answer of it yet after research. And the code isn’t working with pyMC.
I have a forward code where initially I have two known array a and b. I wish to solve an inverse problem and find a probable distribution for a_s and b_s. I am attaching a short layout of my code and my try to implement pyMC in the code.
In my code, the only function that depends on a_s, b_s and time t is the inflow function. With the time integration, q updates. Finally the generated data “bf_rate” is an array. Now using this array data, i want to inversly solve the initial a and b array.
#-------------------------------------------------------------------------------!
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import arviz as az
#-------------------------------------------------------------------------------!
# Input data
#-------------------------------------------------------------------------------!
def Input(--something--):
return -------
#-----------------------------------------------------------------------------!
# Numerical scheme
#-----------------------------------------------------------------------------!
def LaxWendroff(------):
return q
#------------------------------------------------------------------------!
# Given data function
#------------------------------------------------------------------------!
def inflow(----a,b,t,----):
r = 0.5*a[0]
for m in range(1,11):
r = r + a[m]*(np.cos(2.0*(float(m)*np.pi*t/tp))) + b[m]*(np.sin(2.0*(float(m)*np.pi*t/tp)))
r = r*1.0e-6
q[0,j,i] = q[0,j,i]+dt/h[j]*(rate-q[1,j,i+1]*q[0,j,i+1])
q[1,j,i] = rate/q[0,j,i]
return q
#------------------------------------------------------------------------!
# Terminal
#------------------------------------------------------------------------!
def terminal(-------):
return q
#------------------------------------------------------------------------!
# Joint
#------------------------------------------------------------------------!
def joint(-------):
return q
#---------------------------------------------------------------------!
# Output
#---------------------------------------------------------------------!
def output(--------):
return bf
#---------------------------------------------------------------------!
# Defining model for pyMC
#---------------------------------------------------------------------!
def model(a,b):
N,h,A0,c0,CON,rs,rc,iT,CFL = Input(nq,nj,dx)
#---------------------------------------------!
# Time Integration
#---------------------------------------------!
t = 0.0
nt = 10000
nc =7
ntot = nt*nc
t_all = [t]
q_all = [q]
bf_all= [bf]
for k in range(ntot):
t = t + dt
LaxWendroff(ne,nq,NPmax,dt,N,h,A0,c0,q)
inflow(ne,nq,NPmax,a,b,t,h,q)
terminal(ne,nq,NPmax,N,iT,h,A0,c0,rs,q,qold)
joint(ne,nq,nj,NPmax,CON,N,A0,c0,q)
output(ne,nq,Torr,p0,A0,c0,N,q,t)
t_all.append(np.copy(t))
q_all.append(np.copy(q))
bf_all.append(np.copy(bf))
t_all = np.array(t_all)
q_all = np.array(q_all)
bf_all= np.array(bf_all)
q_sample = q_all[60000:70001]
t_sample = t_all[60000:70001]
isegment = 51
ival = int(N[isegment]/4) * np.arange(5)
for i in ival:
bf_rate = q_sample[:,0,isegment,i]*q_sample[:,1,isegment,i]*1.0e6
plt.plot(t_sample,bf_rate,'')
print('Code Ended')
return bf_rate
#%%
a = np.array([172,10,20,-10,5,6,9,-1,40,10,13])
b = np.array([None, 10,20,-10,5,6,9,1,40,-10,13])
# data generation
bf_rate_obs = model(a,b) + np.random.normal(0,0.1)
# %%
# define model
basic_model = pm.Model()
with basic_model:
nas = pm.Uniform("nas", -20, 200, shape=11)
nbs = pm.Uniform("nbs", -20, 50, shape=11)
# Expected & Observed values
q = model(nas,nbs)
observed = bf_rate_obs
# Likelihood (sampling distribution) of observations
sigma = pm.HalfNormal("sigma", sd=0.1)
func = pm.Normal("func", mu=q,sd=sigma,observed=observed)
# trace/sample model
with basic_model:
trace = pm.sample(100, chains=1)
#%%
# plot outputs
with basic_model:
pm.plot_trace(trace)
az.plot_posterior(trace,['a_0'])
plt.xticks(fontsize=10)
az.plot_posterior(trace,['a_1'])