For a classical system consisting of N-interacting identical particles in contact with a heat bath, we define the free energy from thermodynamic relations in equilibrium statistical mechanics. Concretely, the temperature dependence of the free energy is determined from the Gibbs-Helmholtz relation, and its volume dependence is determined from the condition that the quasi-static work in a volume change is equal to the free energy change. Now, we argue the free energy difference in a quasi-static decomposition of small thermodynamic systems. We can then determine the N dependence of the free energy, which includes the Gibbs factorial N! in addition to the phase space integration of the Gibbs-Boltzmann factor.