This paper presents a general framework for smoothing parameter estimation in models with regular likelihoods that depend on unknown smooth functions of covariates. The method is numerically stable and convergent, and it allows for quantifying smoothing parameter uncertainty, which addresses a known issue with the Akaike Information Criterion (AIC) for such models. The framework covers a wide range of models, including generalized additive models (GAMs) with non-exponential family responses, models for location-scale and shape, Cox proportional hazards models, and multivariate additive models. The smooth functions are represented using reduced rank spline-like smoothers with quadratic penalties, and model estimation is performed through penalized likelihood maximization. The smoothing parameters are estimated using Laplace approximate marginal likelihood, which involves Newton optimization and implicit differentiation to obtain derivatives of the model coefficients with respect to the smoothing parameters. The paper also introduces a corrected AIC that accounts for smoothing parameter uncertainty, facilitating model selection. The methods are illustrated through various examples and simulation studies, demonstrating their effectiveness and robustness.This paper presents a general framework for smoothing parameter estimation in models with regular likelihoods that depend on unknown smooth functions of covariates. The method is numerically stable and convergent, and it allows for quantifying smoothing parameter uncertainty, which addresses a known issue with the Akaike Information Criterion (AIC) for such models. The framework covers a wide range of models, including generalized additive models (GAMs) with non-exponential family responses, models for location-scale and shape, Cox proportional hazards models, and multivariate additive models. The smooth functions are represented using reduced rank spline-like smoothers with quadratic penalties, and model estimation is performed through penalized likelihood maximization. The smoothing parameters are estimated using Laplace approximate marginal likelihood, which involves Newton optimization and implicit differentiation to obtain derivatives of the model coefficients with respect to the smoothing parameters. The paper also introduces a corrected AIC that accounts for smoothing parameter uncertainty, facilitating model selection. The methods are illustrated through various examples and simulation studies, demonstrating their effectiveness and robustness.