I need to sample from the following function:

p = p1* p2* p3*p4

Where p1 and p2 are Gaussian functions with the same sigma but different means. But both p3 and p4 are equal to 1/(2*(pi)), both constant functions.

p3 and p4 are defined in the interval -pi to pi.

In the output I need the four variables one from each p1,p2,p3 and p4 composing p. But how can I sample the variables of p3 and p4 if the functions are constant?

I been able to sample from p = p1 * p2 with the emcee Python package in the code below.

from random import sample


import numpy as np
from ROOT import *
import emcee
from math import log10

import ROOT as root


def log_mmc(x,etx,ety):
    return log10((1/(4*3.1415**2))*root.TMath.Gaus(x[0],etx,5)*root.TMath.Gaus(x[1],ety,5))
     


ndim, nwalkers = 4, 100
p0 = np.random.randn(nwalkers, ndim)

sampler = emcee.EnsembleSampler(nwalkers, ndim, log_mmc, args=[12,8])
sampler.run_mcmc(p0, 50000)
 

samples = sampler.get_chain(flat=True`)

So p is as defined above in the code, I didn't bother to normalize it because mcmc does not require it so. But I can't find a way to sample the variable x[2] and x[3] from it, and to set bounds in their values which should be -pi to pi.

How can I achieve this?