Uncertainty quantification for large-scale inverse problems remains a challenging task. For linear inverse problems with additive Gaussian noise and Gaussian priors, the posterior is Gaussian but sampling can be challenging, especially for problems with a very large number of unknown parameters (e.g., dynamic inverse problems) and for problems where computation of the square root and inverse of the prior covariance matrix are not feasible. Moreover, for hierarchical problems where several hyperparameters that define the prior and the noise model must be estimated from the data, the posterior distribution may no longer be Gaussian, even if the forward operator is linear. Performing large-scale uncertainty quantification for these hierarchical settings requires new computational techniques. In this work, we consider a hierarchical Bayesian framework and exploit generalized Golub-Kahan based methods to efficiently sample from the posterior distribution. We also estimate the hyperparameters using a maximum a posteriori estimate of the marginalized posterior distribution, where a stochastic average approximation of the objective function and a preconditioned Lanczos method are used to compute efficient approximations of the function and gradient evaluations. We demonstrate the performance of our approach on various static and dynamic tomography problems.