An iterative learning control (ILC) algorithm has been proposed for the shaping of the output probability density function (PDF) of non-Gaussian stochastic systems. Based on the B-spline modeling of the output PDF, the control strategy is implemented in three steps at each iterative batch. Firstly, the distance between the output PDF of the system and the desired PDF is used between each batch to adjust the B-spline functions. The parameters of the B-spline based PDF model are then identified accordingly. Finally, the control input is designed by optimising a quadratic integral performance function within each batch