Matrix-based predictions and complex variables

4

Hi Cluhmann, thank you for the warm welcome :slight_smile:

Appreciate the shout on the observation data. Correct in saying the values take a wide range. I have tweaked the zero-mean random noise value n_std to a much lower value. The problems persist in that the traces and distributions don’t make sense, particularly the way mass1 and mass2 seem to be bound in the trace to a narrow range, and the sigma values are way higher than stated in the HalfNormal definition.

matrix_tensor_trace.PNG
matrix_tensor_trace.PNG1596×884 269 KB

The version of Pymc I am using is v5.8.2.

Regarding the point about wrapping into a single model context. I endeavored to run as a single model context (assuming I understood correctly - as that is new to me :thinking:).

After some googling etc, I have the following:

import pytensor
import numpy
import pytensor.tensor as pt
from pytensor import function
import numpy as np
import pymc as pm
import matplotlib.pyplot as plt
import arviz as az
from pymc import pytensorf

def generate_synthetic_data():
    # Generate synthetic data
    f = np.linspace(0.01, 10, 1000)  # frequencies for calculation
    w = 2 * np.pi * f  # angular frequencies
    m1_true = 1
    m2_true = 2
    M_true = np.array([[m1_true, 0], [0, m2_true]])

    # observation data Z = jwm
    Z_true = np.zeros((2, 2, len(f)))
    for i in range(len(f)):
        Z_true[:, :, i] = np.abs(1 * w[i] * M_true)

    # add noise to the observation
    n_std = 0.05
    Z_obs = np.abs(Z_true + np.random.normal(0, n_std, Z_true.shape))

    return w, Z_obs  # Return w and Z_obs for later use

def create_pymc3_model(w, Z_obs):
    # Create a PyMC3 model to sample masses and sigma
    with pm.Model() as model:
        # Priors for unknown model parameters
        mass1 = pm.Uniform("mass1", 0, 5)
        mass2 = pm.Uniform("mass2", 0, 5)
        sigma = pm.HalfNormal("sigma", 1)

        # symbolic variables and matrix creation
        m1 = pt.dscalar('m1')
        m2 = pt.dscalar('m2')
        mass_mat = pt.stacklists([[mass1, 0], [0, mass2]])

        # loop step function over w to calculate/predict Z with priors
        def step(w, mass_mat):
            return np.abs(1 * w * mass_mat)

        # store results using scan instead of for looping
        results, _ = pytensor.scan(fn=step, sequences=[w], outputs_info=None, non_sequences=[mass_mat])

        # Likelihood term
        likelihood = pm.Normal("likelihood", mu=results.flatten(), sigma=sigma, observed=Z_obs.flatten())

        # Run the MCMC sampling
        trace = pm.sample(1000, tune=1000)
        
        return trace 

def main():
    plt.close('all')
    
    #main section
    w, Z_obs = generate_synthetic_data()
    trace = create_pymc3_model(w, Z_obs)
   
    #plotting
    az.plot_trace(trace, combined = True)
    
if __name__ == "__main__":
    main()

This runs without crashing however the returned traces remain as per the figure above. I’d be grateful if you can spot anything or see any other issues.

:slight_smile: mmn.