Approximating confidence regions is a central task in parameter estimation, particularly for models governed by differential equations, where nonlinearities limit the accuracy of classical linearized methods. While computationally efficient, these approaches often fail to capture the true geometry of parameter uncertainty. This work presents an extension of the method introduced in [1] for constructing nonlinear confidence regions by iteratively combining local linearizations with likelihood ratio testing. The approach updates the confidence region along principal directions identified via covariance eigenvalue decomposition and adaptively rescales locally linearized confidence ellipsoids to build a Gaussian Mixture Model representation. This yields a probabilistic approximation of the posterior distribution while significantly reducing computational cost compared to full likelihood contouring or sampling-based methods. Numerical experiments demonstrate that the proposed approach accurately captures nonlinear uncertainty structures and closely matches results obtained from Markov Chain Monte Carlo simulations at a fraction of the computational cost. The resulting confidence regions provide a quantitative measure of parameter uncertainty and can be readily integrated into optimal experimental design to guide the selection of informative experiments.