Why not fit 3 large hierarchical models for all subjects, you have appropriate pooling and regularization of parameters from hyperpriors. It is similar to your solution A but not treating each participant independent.
I guess your solution B and C should work as well as an approximation; if all 3 models are equally likely a prior. Some refactoring of the compare function is needed, but the aim is to compare the expected log pointwise predictive density (elpd). In regular use case when you are comparing 3 models independently for each participant, elpd is the same size as the number of observation of that participant (say n trial). So if you are concatenating all subjects together and treating them as one, elpd of each model is now the same size as n_trial*n_subject. I have this work in progress function that might help you a bit.
I think this is valid to do at least as an approximation, but you really should read https://arxiv.org/abs/1704.02030 and do some math/simulation to make sure.