In this paper, a generalized state-space controller design for the shaping of the output probability density function (PDF) is presented for non-Gaussian dynamical stochastic systems. A Radial basis function (RBF) neural network is used to approximate the output PDF of the system. Such a neural network consists of a number of weights and corresponding basis functions. Using such an approximation, the dynamics of the original stochastic system can be expressed as the dynamics between the control input and the weights of the RBF neural network. The task of output PDF control can therefore be reduced to a RBF weight control together with an adaptive tuning of the basis function parameters (i.e., the centers and widths of the basis functions). To achieve this aim, the control horizon is divided into certain intervals hereinafter called batches. Using these definitions, the whole control strategy consists of three stages, namely (a) Sub-space parameter identification of the dynamic nonlinear model (that relates the control signal to the weights of the RBF neural network); (b) Weight tracking controller design using an LMI-based convex optimization technique; and (c) RBF basis functions shape tuning in terms of their centers and widths using an Iterative Learning Control (ILC) framework. Among the above stages, the first two are performed within each batch, while stage (c) is carried out between any two adjacent batches. Such an algorithm has the advantage of the batch-by-batch improvement of the closed-loop output PDF tracking performance. Moreover, the controller mentioned in stage (b) is a general controller in a state-space form. Stability analysis has been performed and simulation results are included to show the effectiveness of the proposed method, where encouraging results have been made.