Inundation, wave arrival and deformation data gathered from measurement campaigns and historical accounts have been used to study past earthquakes and tsunamis, being crucial for magnitude and extension estimations. Nevertheless, the most commonly used methods cannot provide more information about the spatial characteristics of the event's slip distribution. In this work, we aim to create a methodology for obtaining more realistic, heterogeneous estimations of the slip distribution characteristics of past tsunamigenic earthquakes, by obtaining slip patterns of each earthquake of the same seismic segment, thus establishing a history of the seismic cycle in a rupture zone. To obtain the stochastic characterization of the slip distribution of the seismic source, we combine fault defining parameters using a logic tree structure to generate random slip distributions which are subsequently submitted to successive restrictions to assess their compliance with the available deformation and tsunami data, discarding those that do not satisfy the constraints and using those that do to estimate the most probable seismic source in terms of the data. We test this methodology with synthetic heterogeneous slip fields off the coast of South Central Chile to assess its limitations, data dependence, resolution capabilities and uncertainties. Obtaining a good correlation between the synthetic observations and the locations of the main slip features. Finally, the 9.5 Mw, 1960 Valdivia earthquake is used to benchmark the methodology with real deformation and tsunami data from surveys and historical accounts. This result is then compared to slip distributions obtained with tsunami, seismological and deformation data inversions. This shows that the magnitude of the earthquake could be recovered correctly and the slip distribution is very similar to models calculated from different techniques.