Classical methods to infer focal mechanisms from first P wave motions are of limited use for the analysis of small magnitude events occurring in volcanic areas, due to structure complexity. This work presents an alternative way to compute focal mechanisms of earthquakes in volcanic areas, based on a Bayesian approach. The probability density on the whole space spanned by source parameters is graphically represented, by a method originally due to Bernard and Zollo (1989). The data used to infer the probability density function can be of various kinds, like P wave polarities, S wave polarizations and S/P direct wave amplitude ratios. S wave polarizations are left practically unchanged by smooth variations of velocity, and are only sensitive to strong discontinuities and/or seismic anisotropy. S/P amplitude ratios of direct waves are also good observables, not very sensitive to frequency independent amplification effects. They are still affected by differential S-P anelastic attenuation, although to a much lesser extent than absolute S and P amplitudes. In fact, the use of complete waveforms is not desiderable in volcanic areas because wave amplitudes may be strongly biased by not well known path and site effects, mainly for non direct waves which are most affected by medium heterogeneities. Such unmodeled elastic and anelastic features of the medium can be traded off for source characteristics, leading to bad results. The use of the probabilistic approach, and the seloction of robust observables, allow us to solve many problems of classical methods. The determination of the whole probability density on the parameter space visualizes all the information contained in the data set, checking completely the goodness of the solution, and its univoqueness. Furthermore, the use of different kinds of robust data sets helps to constrain the solution, and to minimize the influence of unmodeled medium heterogeneities. In the present method, it is possible to parameterize the source in terms of both double couple sources, and other kinds of sources generally hypothesized on volcanoes (tensile crack, CLVD, explosion) or generic sources characterized by five normalized independent components of the moment tensor. S/P amplitude ratios are shown to be those which most constrain the source mechanism. Furthermore, this paper analyzes the influence of anelastic attenuation on S and P absolute amplitudes, and on S/P amplitude ratios. An example of application to a small (ML = 2.0) earthquake which occurred at Campi Flegrei caldera is reported, showing how the different data sets constrain the solution. Finally, an example is reported of extensive application of the method to the analysis of small earthquakes which occurred at Campi Flegrei (Southern Italy) caldera during an unrest episode (1982-1984). It is shown how the accurate analysis of small earthquakes within a caldera has helped to solve important volcanological questions like: the caldera structure, the mechanism of earthquake generation during unrest episodes, and the features of ground deformations in calderas.