Feature Request: Support for `units` objects in `rvar` data type

I often use the units package to keep track of measurement units and I think it may be interesting to support units objects with the rvar class. I wanted to know whether there is interest to include this feature into the posterior package.

Currently, it is possible to create rvar objects where the draws attribute is a units object:

suppressPackageStartupMessages(library(posterior))
#> Warning: package 'posterior' was built under R version 4.3.2
suppressPackageStartupMessages(library(units))

a <- rvar(units::set_units(1:5, "g"), dim = 1)
draws_of(a)
#> Units: [g]
#>   [,1]
#> 1    1
#> 2    2
#> 3    3
#> 4    4
#> 5    5
class(draws_of(a))
#> [1] "units"

^{Created on 2024-12-10 with reprex v2.0.2}

And simple arithmetic operations work:

class(draws_of(a + a))
#> [1] "units"

But some operations drop the units class. For example:

class(mean(a)) # should return a units object, but returns a numeric vector
#> [1] "numeric"
class(draws_of(rvar_mean(a))) # should keep the units class, but drops it
#> [1] "matrix" "array"

I only did a couple of quick tests, so I have no complete overview, which operations are supported and which are not supported. I can also imagine that adapting some functions may be more difficult.

Start by reviewing the rvar operations mentioned in the report, especially mean(), rvar_mean(), draws_of(), and arithmetic with units objects. Determine which operations should preserve the units class and define the supported behavior; done means the agreed operations retain compatible units and are covered by tests.