A probabilistic approach is presented for jointly inverting gravity gradient and magnetic data for 3-D subsurface distributions of density and magnetic susceptibility. The coupling of the physical property models is incorporated in the inversion by using a cross-covariance matrix of density and magnetic susceptibility. This enables structural similarity such as the orientation and spatial extent of the sources and cross-variance between the two physical properties to be incorporated. A user-defined correlation coefficient can control the level of similarity between the two models. By applying a marginalizing algorithm in the joint inversion, the inversion domain is allowed to be partitioned into various zones, each of which can have its own covariance, cross-covariance matrix, as well as correlation coefficient, depending upon the feature and similarity of sources. Thus, sources with different shapes, sizes and relationships between the two physical properties can be simultaneously recovered. The validity of the method is verified using three synthetic examples, which demonstrate how incorrect parameters of the cross-covariance matrix affect the inverted results. Finally, the proposed method is successfully applied to full tensor gradiometry and magnetic data collected over the Budgell Harbour Stock (BHS) intrusion in north-central Newfoundland, Canada. Compared with models generated from independent inversions, better definition and localization of the main intrusion, as well as associated lamprophyre dykes at shallow depth, are achieved by using the joint inversion. The resolved physical properties for the intrusions show good agreement with field observations of lamprophyre dykes in proximity to the BHS.