Thank your suggestion, Martin! Let me add more background introduction.

The question arised from our recently completed research work implemented by R. When I tried to practise in Stan, it becomes problematic. For the heteroscedastic regression model with unspecified error distribution F_\epsilon that works on aberrant data,

y = \beta^T x + \exp(\gamma^T x_i)\epsilon,

we aim to trim both sides. Let 0 < \tau_1<1/2<\tau_2<1 and F_{\epsilon}^{-1}(\tau) be the lower \tau th quantile of F_{\epsilon}. Then the trimmed mean is written as:

T_{\tau_1, \tau_2}(F_{\epsilon}) = \frac{1}{\tau_2-\tau_1}\int_{F^{-1}(\tau_1)}^{F^{-1}(\tau_2)} sdF_{\epsilon}(s).

We propose the DPM as the prior for f_\epsilon s.t

f_\epsilon (s) =\sum_{k=1}^{\infty} p_k f_{\epsilon_k}(s),

where p_k is drawn from a stick-breaking process and the self-defined kernel is of the following form:

f_{\epsilon_k}(s) = \sum\limits_{l=1}^4 \alpha_{lk}\phi(\frac{s-\mu_{lk}}{\sigma_{lk}}).

Thus it is subject to the constraints and \Phi(\cdot) is the cdf of N(0,1).

Then \alpha is fully determined by (\tau_1, \tau_2) and (\mu_1, \dots, \mu_4, \sigma_1, \dots, \sigma_4).

We also have the prior setting like \mu_{jk} \sim N(0, \sigma_0^2), \sigma_{jk}^2 \sim \text{Inv_Gamma}(\alpha, \beta), j = 1,2,3, 4, k \ge 1.

Notice that \alpha_k should be a weight s.t \sum_{l=1}^{4} \alpha_{lk} =1 and 0 < \alpha_{\cdot k} <1 for any k. That comes the problem that an \alpha determined by the other parameters may not be a legal weight at all! In this case, a resampling scheme is needed to resample the other parameters until the determined \alpha is a legal weight and in R, we use the code we showed in previous question to overcome this problem and that’s why we have to make the posterior not completely deterministic with randomness. But random number generator is not allowed in Stan model block and the resampling scheme is illegal in Stan, which makes us fail to practice that methodology in Stan.

May you give any suggestions to us?