One resource you may like is this recent paper by Culbertson and Sturtz on A Categorical Foundation for Bayesian Probability.
Here are some thoughts on the category $\mathrm{Meas}$ of measurable spaces, where the objects are sets equipped with $\sigma$-algebras, and morphisms are measurable functions (i.e., random variables). This seems like the most natural place to play with foundations of probability, since the category $\mathrm{Meas}$ has a natural concept of a tensor product, where the corresponding $\sigma$-algebra is generated by product sets.
- Question: Is $\mathrm{Meas}$ a Cartesian closed category under the tensor product?
A measure is a countably-additive function defined on $\sigma$-algebras. Typically we learn that measures take values in the real numbers (a Borel-measurable space), but there is no reason that they cannot take values in more general structures, e.g., vector spaces or topological groups. To the best of my knowledge, the most general formulation is due to Tarski, whose monoids are valued in commutative monoids. To Tarski, a "measure" is really a functor from a subcategory of $\mathrm{Meas}$ to the category $\mathrm{ComMon}$ of commutative monoids. This is some universal object, and all other measures can be derived from it. I'm still iffy on this approach, but understanding the Tarski functor is the current thesis project of Tyler Bryson, a masters student at the Courant Institute. He should have more details in a few months.
To define independence, we need to multiply probabilities, so measures should take values in a ring. This multiplicative structure is also important for conditional probabilities, more generally, as well as integration theories. If the space admits symmetries (i.e., a group action), then Tarski's approach spits out a natural multiplicative structure, so I think we're safe in general.
Next,, you can normalize all the measures to have size $1$, turning them into probability measures, but I recommend against it from the categorical point of view. In statistical physics (and more generally, statistics), these normalization constants are the hardest things to compute. I recommend keeping track of them combinatorially, then reducing them at the end. Moreover, it may be useful to keep track of a few different "scales" of measurements, where it is not clear which one to normalize to $1$. This is even seen in the case of vector-valued measures. In the quantum setting of nonnegative-operator-valued measures, a uniform scaling can be chosen, but again the normalization constants are difficult spectral integrals to compute.
No big deal if you want probability measures in the end, just make sure that $1$ makes sense to use in the context you're studying.
Anyway, let $R$ be a topological commutative ring and let $M(X) := M(X,R)$ denote the space of $R$-valued measures on $X$. This is a measurable space, when equipped with the minimal $\sigma$-algebra so that the evaluation map is measurable. This $M$ is an endofunctor on $\mathrm{Meas}$, and closely related to the Giry monad.
Consider the tensor product $XY := X \otimes Y$ of measurable spaces. There are natural projection maps $\pi_X : XY \to X$ and $\pi_Y : X \otimes Y \to Y$. Measures push-forward, so there are natural maps $(\pi_X)_* : M(XY) \to M(X)$ and $(\pi_Y)_* : M(XY) \to M(Y)$. Note that $M(XY)$ corresponds to the joint distributions over $X$ and $Y$, and the resulting projections the marginal distributions.
Note that $M(XY)$ corresponds to the joint distributions over $X$ and $Y$, and the resulting projections the marginal distributions.
There is also the natural space $M(X)M(Y) := M(X) \otimes M(Y)$, and its corresponding projection maps $\pi_{MX} : M(X)M(Y) \to M(X)$ and $\pi_{MY} : M(X)M(Y) \to M(Y)$.
As long as there is some Fubini-type theorem present, there should be a natural map $\varphi : M(X)M(Y) \to M(XY)$ sending a pair of measures to their product measure, corresponding to independent random variables. This should correspond to a commutative diagram, but I don't see it.
Again, watch for the normalization constants, since there are two $M$'s in the source of this map but only one $M$ in the target. Dependent random variables then should be quantified by a failure of some diagram to commute.
Even though I didn't answer the question, hopefully you will find some of this content useful.
If you can improve this community-wiki post, please go for it.