Thanks, this is helpful! I didn’t realize that negloglike = 0.0 would not trigger the same issue for autodiff as the one triggered by pre-allocating using prob = zeros(size(B,1)).
My understanding of the tolerance problem from reading is that the adaptive integtion subdivides the space to do approximate integration until it is considered “good enough” and the number of subdivision depends on my tolerance level. If I switch to the gaussian quadrature, I impose a fixed number of subdivision and the number can be based on a tiny rtol. Then I apply this number to each integration in the loop. It this right?