Thanks! Unfortunately the model is quite complex and pm.find_MAP doesn’t work when I replace the priors with flat priors, however it samples just fine. So can I just record the logp of all the samples and pick the parameter values of the sample with the maximum? Or is there some reason this won’t be the true MLE?