Question

[yebai]

Hi, @nsiccha, README contains the following example

julia_implementation(::Val{:earn_height}; N, earn, height, kwargs...) = begin 
    @stan begin 
        @parameters begin
            beta::vector[2]
            sigma::real(lower=0.)
        end
        @model begin
            earn ~ normal(@broadcasted(beta[1] + beta[2] * height), sigma);
        end
    end
end

Is there a reason for @broadcasted instead of more Julian broadcasting syntax

   earn ~ normal((beta[1] + beta[2] .* height), sigma); 

[nsiccha]

@Broadcasted replaces Base's @., but doesn't materialize the result, ie doesn't allocate. Instead, the broadcasted expression gets fed into the lpdf function (compare

StanBlocks.jl/src/functions.jl

Line 86 in cf7cc37

@inline normal_lpdf(x, mu, sigma) = -bsum(@broadcasted(log(sigma)+.5*square((x-mu)/sigma)))

), where it then gets fed into another (allocation free) @Broadcasted, which then finally gets reduced in a (again allocation free) way that tries to be clever about not doing potentially expensive computation repeatedly (compare the functions in the block starting here:

StanBlocks.jl/src/functions.jl

Line 11 in cf7cc37

begin

).

I think initial Benchmarking revealed that to be faster than doing it the "standard" way. I think that Stan either already does such an optimization OR is conceivably able to do it via Eigen (C++) expressions, ie without the "user" having to do additional stuff like adding the @Broadcasted macro.

Doing it this way, including implementing the lpdf functions myself instead of relying on distributions.jl, made that optimizations easier to do.

[yebai]

Thanks, very helpful.