Uncertainty quantification in geophysical inversion is a well-recognized area of research, yet it has not become routine practice. One of the primary challenges is the computational expense of forward solvers, making robust uncertainty quantification methods like Monte Carlo or Markov Chain Monte Carlo (MCMC) impractical, particularly for higher-dimensional problems. This challenge is amplified in the case of joint inversion, where multiple types of forward solvers must be run thousands of times. We propose a stochastic joint inversion framework that integrates the Very Fast Simulated Annealing (VFSA) approach with a generalized fuzzy c-means clustering technique for effective parameter coupling. By incorporating a sparse parameterization strategy and executing multiple VFSA chains with varying initial models, we effectively mitigate VFSA's tendency to converge at the peak of the derived posterior probability density (PPD) function. The approach presented here address the inherent challenge of high computational costs for implementing joint-inversion with nonlinear sampling methods like MCMC by providing a feasible probabilistic joint inversion alternative that can integrate petrophysical information as well as geological constraints.