A method for parametric estimation of seismic wavelets from well logs and seismic data is developed. Parameters include amplitude, skewness, length and fluctuation order, and the link between parameters and wavelet properties provides a user-friendly interpretation of the wavelet function. The method is set in a Bayesian framework, and is well-suited for addressing questions about uncertainty related to estimated wavelets. This is accomplished by sampling the posterior distribution using Markov Chain Monte Carlo methods. The estimation method is framed as a practical step-wise procedure. An extension of the model to enable joint wavelet estimation from seismic data with multiple incidence angles, is also described. The method is tested on simulated data, and on well log and seismic amplitude data from the North Sea. The results in the synthetic case indicate that the method performs well under idealised conditions. When tested on real data, the method produces a realistic wavelet fit and uncertainty range. Uncertainty is substantially reduced from the prior to the posterior distribution, but in general, the shape of the posterior surface could make it hard to explore. A comparison with a wavelet estimator based on a Gaussian process indicates that the proposed parametric form gives a tighter wavelet, and is less prone to overfitting.