I'd like to understand how functions like Mean are overloaded to provide the right behavior for the distributions in Mathematica (like NormalDistribution or PoissonDistribution).

I originally assumed it was through the use of UpValues, but now I'm not so sure...

So, if I wanted to implement a function or distribution and define a behavior for it when another function like Mean is applied to it, how would I go about it? I know it's not the following:

f[x_]:=2+x
f /: Mean[f[x_]] := 3 x

Desired behavior:

f[3]
Mean[f[3]]

(*
==> 5
==> 9
*)