Problem with using pm.Deterministic to return a matrix when using a function as input to pm.deterministic and having prior as a argument to that function

,
1 
import pymc as pm
import numpy as np
import pandas as pd
import arviz as az
import pytensor.tensor as pt
import scipy.special
from scipy.stats import norm
from scipy.special import logit
import os
import graphviz
import scipy
import time


def create_gamma(T, r, q):
    vector = np.zeros(T)
    vector[0] = (1 - q)
    for k in range(1, T):
        vector[k] = scipy.special.comb(k + r - 1, r - 1) * (1 - q)**(k+1)
    return vector

def create_time_lagged_N(data, n):
    matrix = np.zeros((n, n))  # Ensure matrix has correct dtype
    for i in range(n):
        matrix[i, :i+1] = data[i::-1]
    return matrix

def create_covariance_matrix(T, q):
    # Initialize an empty matrix of size T x T
    matrix = np.zeros((T, T))

    # Define the elements based on the provided pattern
    for i in range(T):
        for j in range(T):
            if i == j:
                # Diagonal elements
                matrix[i, j] = 1 + 4*(1-q)**2 + (1-q)**4
            elif j == i + 1 or i == j + 1:
                matrix[i, j] = -2*(1-q)*(1 + (1-q)**2)
            elif j == i + 2 or i == j + 2:
                matrix[i, j] = (1-q)**2

    return matrix

def lag_function(x, T, q, r):
    N_gamma = create_time_lagged_N(x, T)
    gamma = create_gamma(T, r, q)
    return np.dot(N_gamma, gamma)

def do_MCMC(x, y, T):
    # PyMC model
    with pm.Model() as model:
        # Priors for the parameters
        p = pm.Uniform('p', lower=0, upper=1)
        q = pm.Uniform('q', lower=0, upper=1)
        r = pm.DiscreteUniform('r', lower=1, upper=6)
        sigma2 = pm.InverseGamma('sigma2', alpha=2.0, beta=1.0)

        # Create deterministic covariance matrix V
        V = pm.Deterministic('V', create_covariance_matrix(T, q))

        # Create deterministic lag function
        mean_mret_t = pm.Deterministic('mean_mret_t', p * q ** r * lag_function(x, T, q, r))

        # Likelihood
        likelihood = pm.MvNormal('mret_t', mu=mean_mret_t, cov=sigma2 * V, observed=y)

        # Sampling using Metropolis-Hastings
        step = pm.Metropolis()
        trace = pm.sample(2000, step=step, tune=1000, return_inferencedata=False)

        return trace, model

I am confused about the

V = pm.Deterministic('V', create_covariance_matrix(T, q))

and

mean_mret_t = pm.Deterministic('mean_mret_t', p * q ** r * lag_function(x, T, q, r))

Both are giving errors, I believe something is going wrong when create_covariance_matrix or lag_function gets called with the prior parameters q and r. I have never used such a complicated function such as create_covariance_matrix. Could someone direct me in the right direction?

The error that I am getting is

Exception has occurred: ValueError
setting an array element with a sequence.

for the lines

vector[0] = (1 - q)

and

matrix[i, j] = 1 + 4*(1-q)**2 + (1-q)**4
2

Welcome!

You seem to be mixing numpy arrays and PyMC/PyTensor tensors (and scipy functions). q is a tensor, but vector and matrix are numpy arrays. Numpy doesn’t know what to do with a tensor. You can create empty data structures in PyTensor and things should work mostly the same as with numpy. Check out the getting started and basic tensor functionality pages for relevant documentation. For example, vector should probably be a pt.vector and matrix should probably be a pt.matrix.

1 Like
3

For assignment you’ll have to use pt.set_subtensor.

x[0] = 1 becomes x = pt.set_subtensor(x[0], 1) (you have to use the returned object)

1 Like