Piecewise linear quadratic (PLQ) penalties play a crucial role in many applications, including machine learning, robust statistical inference, sparsity promotion, and inverse problems such as Kalman smoothing. Well known examples of PLQ penalties include the l2, Huber, l1 and Vapnik losses. This paper builds on a dual representation for PLQ penalties known from convex analysis. We provide conditions that allow these losses to be interpreted as negative logs of true probability densities, and enable construction of non-smooth multivariate distributions with specified means and variances from simple scalar building blocks. The main contribution of this paper is a flexible statistical modelling framework for a variety of learning applications, where solutions to all models in this framework can be computed using interior point (IP) methods. IP methods solve nonsmooth optimization problems by working directly with smooth systems of equations characterizing the optimality of these problems. The efficiency of the IP approach depends on the structure of particular applications, and we consider the class of dynamic inverse problems using Kalman smoothing. This class comprises a wide variety of applications, where the aim is to reconstruct the state of a dynamical system with known process and measurement models starting from noisy output samples. In the classical case, Gaussian errors are assumed both in the process and measurement models for such problems. The extended framework allows arbitrary PLQ densities to be used, and the proposed IP approach solves the generalized Kalman smoothing problem while maintaining the linear complexity in the size of the time series, just as in the Gaussian case. This extends the computational efficiency of the Mayne-Fraser and Rauch-Tung-Striebel algorithms to a much broader nonsmooth setting, and includes many recently proposed robust and sparse smoothers as special cases.