This work considers the estimation of impulsive time series pertaining to biomedical systems and, in particular, to endocrine ones. We assume a signal model in the form of the output of a continuous linear time-invariant system driven by a sequence of instantaneous impulses, which concept is utilized here, in particular, for modeling of the male reproductive hormone axis. An estimation method to identify the impulsive sequence and the continuous system dynamics from sampled measurements of the output is proposed. Hinging on thorough mathematical analysis, the method improves upon a previously developed least-squares algorithm by resolving the trade-off between model fit and input sparsity, thus removing the need for manual tuning of user-defined estimation algorithm parameters. Experiments with synthetic data and Markov chain Monte-Carlo estimation demonstrate the viability of the proposed method, but also indicate that measurement noise renders the estimation problem ill-posed, as multiple estimates along a curve in the parameter space yield similar fits to data. The method is furthermore applied to clinical luteinizing hormone data collected from healthy males and, for comparability, one female, with similar results. Comparison between the estimated and theoretical elimination rates, as well as simulation of the estimated models, demonstrate the efficacy of the method. The sensitivity of the impulse distribution to the estimated elimination rates is investigated on a subject-specific data subset, revealing that the input sequence and elimination rate estimates can be interdependent. The dose-dependent effect of a selective gonadotropin releasing hormone receptor antagonist on the frequency and weights of the estimated impulses is also analyzed; a significant impact of the medication on the impulse weights is confirmed. To demonstrate the feasibility of the estimation approach for other hormones with pulsatile secretion, the modeling of cortisol data sets collected from three female adolescents was performed.