How to use numpy array as input in the PyMC model?

Hello everyone,

I am new in PyMC. I am writing a Fourier function to generate some data which has two array of coefficients. Using PyMC, I want to predict the known coefficient at the end. Something like below:

import numpy as np
import matplotlib.pyplot as plt
import arviz as az
import pymc3 as pm
import theano.tensor as tt

def model(a,b,t):
    for i in range(2):
        f = a[i]*(np.cos(2.0*np.pi*t)) + b[i]*(np.sin(2.0*np.pi*t))
    return f

# generate data
n=20
sig = 0.2
t = np.linspace(0,1,n)

# True Values with noise
a = np.array([72,88])          
b = np.array([13,10])
for i in range(2):
    f_obs = a[i]*(np.cos(2.0*np.pi*t)) + b[i]*(np.sin(2.0*np.pi*t)) + np.random.normal(0,0.1)

# define model
basic_model = pm.Model()
with basic_model:
    a_array = np.zeros(2)
    b_array = np.zeros(2)
    a_now = pm.Uniform('a_%d',lower=70, upper=90)
    b_now = pm.Uniform('b_%d',lower=8, upper=15)
    a_array.append(a_now)
    b_array.append(b_now)
    
    a = tt.as_tensor_variable(a_array)
    b = tt.as_tensor_variable(b_array)

    # Expected & Observed values
    q = model(a,b,t)
    observed = f_obs

    # Likelihood (sampling distribution) of observations
    sigma = pm.HalfNormal("sigma", sd=0.01)
    func = pm.Normal("func", mu=q,sd=sigma,observed=observed)

    # trace/sample model
    with basic_model:
        trace = pm.sample(1000, chains=1)

    #%%
    # plot outputs
    with basic_model:
        pm.plot_trace(trace) # pm.plot_trace(trace,['nu','sigma'])
        az.plot_posterior(trace,['a_1'],ref_val=2)
        plt.xticks(fontsize=10)
        az.plot_posterior(trace,['b_1'],ref_val=2)
        plt.xticks(fontsize=10)
        az.plot_posterior(trace,['sigma'],ref_val=np.pi/2)
        #plt.xticks(fontsize=10)

If I take single a and b, it works, But it is not working for numpy array of a and b. I am not sure about the PyMC part. Can anyone direct me to my errors?

First, I’m not an expert, but…
Your f and fobs look like they would only be based on index of value = 1.
Inside of model() function it seems like you would want your entire tensor operated on (remove the for loop and indexing)
Try using graphviz to look at the dimensionality of your observed and your a,b,q elements.

Can you share the error messages?