In this paper, an Iterative Learning Control (ILC) scheme is presented for the control of the shape of the output probability density functions (PDF) for a class of stochastic systems in which the relationship between approximation basis functions and the control input is linear, and the stochastic system is not necessarily Gaussian. A Radial Basis Function Neural Network (RBFNN) has been employed for the output PDF approximation and the coefficients of the approximation are linearly related to the control input. A three-stage method for the ILC-based PDF control is proposed which incorporates a) identifying PDF model parameters; b) calculating the control input; and c) updating RFBN parameters. The latter is accomplished based on P-type ILC law and the difference of the desired and calculated output PDF within a batch. Conditions for the convergent ILC rules have been derived. Simulation results are included to demonstrate the effectiveness of proposed method.