This is true. I actually realize now that this example of a ‘complicated distribution’ was probably not the best example since it is stated on the page where I took the code (see first post) that “The data only enters via […] which is completely independent of the fit parameter alpha and can be computed once at beginning. Therefore, handling huge datasets such as a Gaia catalogue is no problem at all!”
This is definitely not the path I took, since as you mention, the function sampleFromSalpeterDATA is called each time, hence the time needed to compute. In the original code they use the logp, which in that case does not depend on the fit parameter alpha.
This brings me to ask, like I mentioned in the first post, if there are examples out there of people running external codes with PyMC, just using Metropolis (since no gradients are then available). I just found this and this so far.
Is it even the best idea to use PyMC to just use Metropolis-Hastings on an external model?