I dont think I know (in my limited experience) how to do it starting from the model that you already have but if I understand your question I would simply build Pmarginal myself as histogram and then sample data from there. Never used Ricardo’s suggestion the do calculus though that also seems like an interesting option. In more detail:
1- Since you know the exact forms of P1,P2 and your priors are uniform (they are basically just setting boundaries for your integrals), I would check if the integral you have given can be calculated exactly (if you are lucky you might get something in terms of incomplete gamma functions). If not I would integrate it numerically. In both cases you get a density (up to scale) for Pmarginal depending on ng. Using that density as an histogram you can generate random Ng observables, see for instance:
2- Then you use these randomly generated observables just as you have used ng in your model. This being said, I have never attempted such a thing myself so I don’t have any feeling on how the resolution you choose for your integration and histogram sampling would affect the results you would get, so that should be something that is also tested.
ps:
Note that in pymc_experimental there is marginalization over finite discrete parameters
https://www.pymc.io/projects/experimental/en/latest/generated/pymc_experimental.MarginalModel.html
but I suspect that would not be sufficient for you since you want to marginalize over continuous parameters. But still putting it here in case you find it useful somehow.