What exactly are the characteristics of the built-in, vectorized functions that make them superior to coding log-likelihoods and adding to target by hand? Is it that the analytic derivatives are coded up behind the scenes, making the differentiation faster? superior numeric stability? prevention of loops via matrix math? It’s probably a combination of all three, but how should I think about what facets of my particular problem will influence the relative importance of these considerations?
For background, I am working with a likelihood involving censored Weibull-like variables. I do not believe this likelihood can be vectorized with weibull_lpdf and weibull_lccdf because those functions return the log-pdf and log-survival summed across the entire input vector. (I know I haven’t explained well why I can’t vectorize, but I don’t want to get too bogged down in the specifics of my particular problem. If you are unusually curious or just like reading ugly code, the Stan file is here.)
Right now my code loops through the N observations, calling the weibull_lpdf and weibull_lccdf functions when possible in order to obtain component likelihood quantities like the log-instantaneous hazard. Are there ways to do DIY vectorization to get most of the speed/stability benefits of the built-in functions? Should I be writing my own _lpdf function and somehow specifying analytic derivatives somewhere?
You can’t do it from within the Stan language. You have to declare but not define a user-defined function and #include a C++ definition for that function. But in the case of just looping over a vector to get component-wise likelihoods, I think calling weibull_lpdf and weibull_lccdf is fine.
See the discussion in the manual chapter on efficiency for the full story.
The biggest issue is a compressed autodiff graph leading to fewer virtual function calls and more memory locality in the backward pass of autodiff.
The second big issue is shared computations, i.e., only computing log(sigma) once in y ~ normal(mu, sigma) when y is a vector.
The third big issue is that we drop constant terms using the ~ notation, so y ~ normal(mu, sigma) is more efficient than target += normal_lpdf(y | mu, sigma).