My model has 3 independent parameters $\{\rho, \alpha, \beta\}$ (polar coordinates), and a set of observables $\{Q_i\}$ and $\{T_{ij}\}$ where $i=1,2,...,642$ and $j=1,2,...,Q_i$ (if $Q_i=0$, there is no $T_{ij}$).
From $\{\rho, \alpha, \beta\}$, a set of intermediate (spherical harmonic) coefficients $q^m_l(\rho, \alpha, \beta)$ and $t^m_1(\rho, \alpha, \beta)$ where $l=1,2,3~m=-l,...,-1,0,1,...,l$ are obtained:
A rotational transformation $R_l^m(\alpha, \beta)$ (trigonometric polynomials in $(\alpha, \beta)$):
$q^m_l(\rho, \alpha, \beta) = q^0_l(\rho, 0, 0) R^m_l(\alpha, \beta)$
$t^m_1(\rho, \alpha, \beta) = t^0_1(\rho, 0, 0) R^m_l(\alpha, \beta)$
$q^0_l(\rho, 0, 0)$ and $t^0_l(\rho, 0, 0)$ are calibrated empirical functions.
$Q_i$ obeys Poisson distribution[a],
$Q_i \sim \text{Pois}(\sum_{l=1}^3\sum_{m=-l}^lq^m_l(\rho, \alpha, \beta)S^m_l(\theta_i, \phi_i))$, where $S^m_l(\theta_i, \phi_i)$ is the value of real spherical harimonic function at $i$.
$T_{ij}$ obeys shifted exponential and is treated with quantile regression[b] at 1% (empirical) percentile,
$T_{ij} \sim \text{Quant}_{0.01}(\sum_{l=1}^3\sum_{m=-l}^lq^m_l(\rho, \alpha, \beta)S^m_l(\theta_i, \phi_i))$
I would like to do inference, and I can see 2 options.
Hierarchical Bayesian model: use the above generative process and fit with MCMC. I know stan and pymc.
I have millions of observations, and am concerned about the feasibility of this approach on a small cluster (less than 1000 processor cores).
Ensemble of predictors: given an observation $\{Q_i, T_{ij}\}$, fit $\hat{q}^m_l$ and $\hat{t}^m_l$. For each $l$ fit $\rho, \alpha, \beta$ by $q$ or $t$. With 6 predictors $(\hat{\rho}_l^{t/q}, \hat{\alpha}_l^{t/q}, \hat{\beta}_l^{t/q})$, what is the best way to make a combined predictor? I know boosting for classifiers, is there a continuous counterpart?
Am I going to lose accuracy compared to hierarchical Bayesian? How can I estimate that loss?
Footnotes:
a. actually quasi-Poisson, where the variance is proportional to mean.
b. Can a generalized linear model use shifted exponential as residual distribution?