This paper presents a finite size scaling analysis of the block distribution functions of the Ising model in dimensions d=2, 3, and 4. The distribution function $ P_L(s) $ of local order parameters in blocks of linear dimension L is studied for both subsystems of infinite lattices and finite systems with free and periodic boundary conditions. Above the critical temperature $ T_c $, these distributions tend toward a Gaussian centered at zero block magnetization, while below $ T_c $, they tend toward two Gaussians centered at $ \pm M $, where M is the spontaneous magnetization. However, the wings of the distribution at small $ |s| $ are non-Gaussian, reflecting two-phase coexistence. These distribution functions can be used to estimate interface tension between ordered phases. At criticality, the distribution functions tend toward scaled universal forms, dependent on boundary conditions. These scaling functions are estimated from Monte Carlo simulations. For subsystem-blocks, good agreement is found with previous renormalization group work. The paper shows that Monte Carlo studies of critical phenomena can be improved using these distribution functions. They can improve estimates of order parameter, susceptibility, and interface tension. $ T_c $ can be estimated independently of critical exponent estimates. A Monte Carlo "renormalization group" similar to Nightingale's phenomenological renormalization is proposed, yielding accurate exponent estimates. Information on coarse-grained Hamiltonians can be gained, which is particularly interesting for more general Hamiltonians. The paper introduces the concept of dividing systems into blocks of finite linear dimension L. For a d-dimensional Ising system, the magnetization $ s_i $ of the i-th block is defined. The Boltzmann factor is replaced by the probability for the field $ s_i $, leading to a Ginzburg-Landau-Wilson Hamiltonian. The parameters $ r_L, u_L, v_L, \ldots, C_L, \ldots $ are related to the microscopic Hamiltonian. However, the resulting theories can only predict universal properties, losing information on non-universal properties. Monte Carlo methods are used to obtain numerical results by sampling the distribution function $ P_L(\{s_i\}) $. The paper studies the reduced distribution function of one block and related distribution functions for finite systems. Section II describes general properties of these functions, including properties for $ L \gg \xi $, a situation not normally considered in the renormalization group approach. Section III presents numerical results and universal scaling functions. Section IV proposes a phenomenological "Monte Carlo renormalization group (MCRG)" similar to Nightingale's finite size renormalization group. Section V contains conclusions and discusses generalizations to other systems.This paper presents a finite size scaling analysis of the block distribution functions of the Ising model in dimensions d=2, 3, and 4. The distribution function $ P_L(s) $ of local order parameters in blocks of linear dimension L is studied for both subsystems of infinite lattices and finite systems with free and periodic boundary conditions. Above the critical temperature $ T_c $, these distributions tend toward a Gaussian centered at zero block magnetization, while below $ T_c $, they tend toward two Gaussians centered at $ \pm M $, where M is the spontaneous magnetization. However, the wings of the distribution at small $ |s| $ are non-Gaussian, reflecting two-phase coexistence. These distribution functions can be used to estimate interface tension between ordered phases. At criticality, the distribution functions tend toward scaled universal forms, dependent on boundary conditions. These scaling functions are estimated from Monte Carlo simulations. For subsystem-blocks, good agreement is found with previous renormalization group work. The paper shows that Monte Carlo studies of critical phenomena can be improved using these distribution functions. They can improve estimates of order parameter, susceptibility, and interface tension. $ T_c $ can be estimated independently of critical exponent estimates. A Monte Carlo "renormalization group" similar to Nightingale's phenomenological renormalization is proposed, yielding accurate exponent estimates. Information on coarse-grained Hamiltonians can be gained, which is particularly interesting for more general Hamiltonians. The paper introduces the concept of dividing systems into blocks of finite linear dimension L. For a d-dimensional Ising system, the magnetization $ s_i $ of the i-th block is defined. The Boltzmann factor is replaced by the probability for the field $ s_i $, leading to a Ginzburg-Landau-Wilson Hamiltonian. The parameters $ r_L, u_L, v_L, \ldots, C_L, \ldots $ are related to the microscopic Hamiltonian. However, the resulting theories can only predict universal properties, losing information on non-universal properties. Monte Carlo methods are used to obtain numerical results by sampling the distribution function $ P_L(\{s_i\}) $. The paper studies the reduced distribution function of one block and related distribution functions for finite systems. Section II describes general properties of these functions, including properties for $ L \gg \xi $, a situation not normally considered in the renormalization group approach. Section III presents numerical results and universal scaling functions. Section IV proposes a phenomenological "Monte Carlo renormalization group (MCRG)" similar to Nightingale's finite size renormalization group. Section V contains conclusions and discusses generalizations to other systems.