Solving Inverse problem using pyMC

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'])

Can anyone help?

Hi, I’m also very new to pyMC, so take my input with a big grain of salt.

I think you may have issues inputing pyMC3 variables (nas,nbs) directly into your function.
In addition, if your model returns a vector, func should perhaps be multivariate.

In all cases, I’m not sure you can avoir using a black-box likelihood here, but I may be wrong. There is an example notebook here : Using a “black box” likelihood function — PyMC3 3.11.5 documentation

I hope this helps

1 Like

Hi @salmanromeo40
So in theory you can just write an arbitrary generative model and use that to sample. In practice, there needs to be a way to calculate gradient information (not just calculate log probabilities). In PyMC models, there this is done by building up a symbolic graph in Aesara and that gives us the gradients. There is some level of numpy/python code that can be included in PyMC models which can be ‘understood’ and used to generate Aesara graphs, but there is a certain point at which this breaks. The exact point, I do not fully understand (maybe @ricardoV94 can shed light on this).

But your model seems to include multiple (partially redacted) functions. If these are complex then it might not be possible for it to be parsed into an Aesara graph and so the gradient information isn’t available. In these cases, the general approach would be to write your functions using Aesara operations, rather than raw Python or numpy operations.

2 Likes