Gibbs' analysis (see his 1902 book, and the beautifully concise summary of his main derivations in Chapter 10 of The Logic of Thermostatistical Physics by Emch and Liu) was purely mathematical based on the assumption of an underlying Hamiltonian description of the system, using the phase space density $\rho(X)$, where $X$ is a point in phase space. Within this analysis he found that the expectation of the log of the phase space density, $\bar{\eta} = \langle \eta \rangle$ — where $\eta = \log(\rho)$, so that:
$$
\langle \eta \rangle= \int \rho(X) \, \log\left[\rho(X)\right] \, {\rm d} \Gamma
$$
with ${\rm d}\Gamma$ an infinitesimal volume element in the phase space — played the same role as (negative) entropy in a derived equation that matched the fundamental thermodynamic relation
$$
{\rm d} U = T \, {\rm d}S - P\, {\rm d}V
$$
He therefore proposed that his completely mechanistic definition be considered an expression for thermodynamic entropy.
See Chapter IV in his book (linked above), Equation (114) and the preceding discussion; and Chapter XIV, surrounding equation (483). (Forewarning: it is difficult to parse because he uses the word "index" for "log" and then uses log[exp(x)] = x to write things in an unfamiliar way).
Using the Gibbs Entropy as a functional over distributions $\rho$, you could determine the canonical distributions by extremizing with constraints. This is accomplished by calculus of variations, but it is the same idea, and easier to see, using Shannon's discrete version:
$$
S = -k \sum_i p_i \log p_i
$$
- If we assume only that $\sum p_i = 1$, then the multi-dimensional extremization of the multivariable function $S\left(\left\{p_i\right\} \right)$, subject to the constraint that $\sum p_i$ is constant, is accomplished with a Lagrange multiplier, $\lambda$:
$$
\frac{\partial}{\partial p_j} \left[ S\left(\left\{p_i\right\} \right) + \lambda \sum_i p_i \right] = 0
$$
which leads to $p_i = 1/N$ for all $i=1,\ldots, N$. This is the microcanonical distribution.
- If we additionally assume a constant expectation of the energy, $U = \sum_i \varepsilon_i p_i$, then:
$$
\frac{\partial}{\partial p_j} \left[ S\left(\left\{p_i\right\} \right) + \lambda_1 \sum_i p_i + \lambda_2 \sum_i \varepsilon_i p_i \right] = 0
$$
leads to $p_i = A \, \exp\left(- B \, \varepsilon_i\right)$, which is the canonical distribution.
- The same procedure can be used to determine the grand canonical distribution, by additionally assuming a fixed expected particle number, $N = \sum n_i p_i$.
Finally (as Emch and Liu make clear in their book), if we consider a system with distribution $\rho_0$ on an aggregate phase space, $\Gamma = \Gamma_1 \times \Gamma_2$, composed of two subsystems, and define the reduced probability densities of each sub-system:
$$
\rho_1 = \int_{\Gamma_2} \rho_0 \quad \text{and} \quad \rho_2 = \int_{\Gamma_1} \rho_0
$$
then the corresponding entropies
$$
S_n = -k \int_{\Gamma_n} \rho_n \log \rho_n \quad (n=0,1,2)
$$
satisfy the inequality:
$$
S_0 \le S_1 + S_2
$$
with
$$
S_0 = S_1 + S_2 \quad \text{iff} \quad \rho_0 = \rho_1 \cdot \rho_2
$$
i.e., if the subsystems are independent.
Jaynes extended these facts to "prove" the Second Law of Thermodynamics. See his paper Gibbs vs Boltzmann entropies (1965) and this physics.se answer.
