Critical finite size scaling functions for the order parameter distribution of the two and three dimensional Ising model are investigated. Within a recently introduced classification theory of phase transitions the universal part of critical finite size scaling functions has been derived by employing a scaling limit which differs from the traditional finite size scaling limit. In this paper the analytical predictions are compared with Monte Carlo simulation results. We find good agreement between the analytical expression and the simulation results. The agreement is consistent with the possibility that the functional form of the critical finite size scaling function for the order parameter distribution is determined uniquely by only a few universal parameters, most notably the equation of state exponent.