Is there a Lagrangian formulation of statistical mechanics? In statistical mechanics, we usually think in terms of the Hamiltonian formalism. At a particular time $t$, the system is in a particular state, where "state" means the generalised coordinates and momenta for a potentially very large number of particles. (I'm interested primarily in classical systems for the sake of this question.) Since this state cannot be known precisely, we consider an ensemble of systems. By integrating each point in this ensemble forward in time (or, more often, by considering what would happen if we were able to perform such an integral), we deduce results about the ensemble's macroscopic behaviour. Using the Hamiltonian formalism is useful in particular because it gives us the concept of phase space volume, which is conserved under time evolution for an isolated system.
It seems to me that we could also consider ensembles within the Lagrangian formalism. In this case we would have a probability distribution over initial values of the coordinates (but not their velocities), and another distribution over the final values of the coordinates (but not their velocities). (Actually I guess these would need to be two jointly distributed random variables, since there could easily be correlations between the two.) This would then lead to a probability distribution over the paths the system takes to get from one to the other. I have never seen this Lagrangian approach mentioned in statistical mechanics. I'm curious about whether the idea has been pursued, and whether it leads to any useful results. In particular, I'm interested in whether the idea of phase space volume has any direct meaning in terms of such a Lagrangian ensemble.
 A: There is a field theory version of statistical physics. The temperature is like the imaginary time. In this way we can formulate theory by path integral with action determined by Lagrangian.
A: I am not sure if this is what you are up to (it is related to what Xiao-Qi Sun said) to but I'll give it a try too ...
At the beginning of Chapter V.2 of his QFT Nutshell, Anthony Zee explains how classical statistical mechanics (characterized by the corresponding partition function involving the Hamilton function) in $d$- dimensional space is related to Eucledian field theory (characterized by the corresponding generating functional or path integral involving the Lagrangian). 
To see this relationship, consider for example the Minkowskian path integral of a scalar field
$$
(1) \,\,
\cal{Z} = \int\cal{D}\phi e^{(i/\hbar)\int d^dx[\frac{1}{2}(\partial\phi)^2-V(\phi)]}
        = \int\cal{D}\phi e^{(i/\hbar)\int d^dx\cal{L}(\phi)} 
        = \int\cal{D}\phi e^{(i/\hbar)S(\phi)}
$$
Upon Wick rotation, the Lagrange density $\cal{L}(\phi)$ turns into the energy density and the action $S(\phi)$ gets replaced by the energy functional $\cal E(\phi)$ of the field $\phi$
$$
(2) \,\,
\cal{Z} = \int\cal{D}\phi e^{(-1/\hbar)\int d^d_Ex[\frac{1}{2}(\partial\phi)^2+V(\phi)]}
        = \int\cal{D}\phi e^{(-1/\hbar)\cal{E}(\phi)} 
$$
with
$$
\cal E(\phi) = \int d^d_Ex[\frac{1}{2}(\partial\phi)^2+V(\phi)]
$$
This can now be compared to the classical statistical mechanics of an N-particle system with the Energy
$$
E(p,q) = \sum_i \frac{1}{2m}p_i^2+V(q_1,q_2,\cdots,q_N)
$$
and the corresponding partition function 
$$
Z = \prod_i\int dp_i dq_i e^{-\beta E(p,q)}
$$
Integrating over the momenta $p_i$ one obtains the reduced partition function 
$$
Z = \prod_i\int dq_i e^{-\beta V(q_1,q_2,\cdots,q_N)}
$$
Following the usual procedure to obtain the field theory which corresponds to this reduced partiction function by letting $i\rightarrow x$, $q_i \rightarrow \phi(x)$ and identifying $\hbar = 1/\beta = k_B T$ it has exactly the same form as the Euclidian path integral (2).
So it can finally be seen that in this example, the (reduced) partition function of an N-particle system in d-dimensional space corresponds to the path integral of a scaler field in d-dimensional spacetime.
These arguments can be further generalized to obtain a path integral representation of the quantum partition funcction, finite temperature Feynman diagrams, etc too ...
If I understand this right, this line of thought relating statistical mechanics to field theory is for example applied in topics like the Nonequilibrium functional renormalization group or in AdS/CFT to relate the correlation functions on the CFT side to the string amplitudes on the AdS side.
A: The Hamilton formulation of classical dynamics gives rise to a very strong and important theorem in statistical mechanics that is the Liouville theorem. As you probably know already it states that the probability density $\rho(\mathbf{r}, \mathbf{p})$ to be around a given point $(\mathbf{r}, \mathbf{p})$ in phase space follows the equation of evolution:
$\frac{\partial \rho}{\partial t} = \{\rho, H \}$ where $\{ \cdot\}$ denotes the Poisson brackets.
This equation is equivalent to Hamilton equations of evolution for $(\mathbf{r}, \mathbf{p})$.
Now, when you look at macrovariables, it can be worked out (it has been done first by Zwandsig I think) that the Liouville equation (for the microvariables) gives rise to a Fokker-Planck equation for these macrovariables. It is in spirit very similar to the Liouville equation except that there is a stochastic component in it whose simplest characteristic is to add a second space derivative on the right hand side of the evolution equation.
Now, if you know your maths, you also know that any Fokker-Planck equation can be associated to a set of stochastic equations for the macrovariables under study (one very famous being the Langevin equation)...and we are back to something very close to the Hamilton equations but for macrovariables.
In case you were wondering if there is a minimum action principle for these stochastic equations, I am not aware of that. I think, they are very similar to Shrodinger equation in this respect. However what it means is that indeed the macrovariable propagators can be expressed as path integrals. The Wiener measure is one typical case.
Note that my answer is focused on Hamilton and Lagrangian dynamics in the classical sense where they were used to compute trajectories in time.
In classical statistical mechanics, you could find a Lagrangian approach akin to what is done in, say, QFT. This would be the Landau-Ginsburg approach of phase transitions and complex systems in general.
A: I'm not sure this is quite what you're looking for, as classical statistical mechanics in equilibrium is, by definition, not a properly dynamical system, but min/maxing the function $$S[p_i] = -k_{B}\sum_{i}{\ln{p(x)}} + Z\sum_i{\big(p(x)-1\big)}+\beta\sum_i{p_iE_i-\langle E\rangle}+\mu\sum_i{N_ip_i-\langle N\rangle}+...$$
not only recovers the Boltzmann distribution, but naturally gives rise to the second law of thermodynamics and provides a unified way of dealing with any kind of ensemble. As such, this S is occasionally described as the Action for an equilibrium distribution
A: The transition between the Hamiltonian and Lagrangian formalisms in Mechanics can be accomplished by means of the Hamilton-Jacobi theory.
Consider for example a classical statistical ensemble on a phase space $(x,p)$ defined by:
A. The (initial) state of this ensemble is defined by a distribution function $f(x_0,p_0)$  satisfying the normalization condition:
$$\displaystyle{\int f(x_0,p_0) dx_0dp_0 = 1}$$
($(x_0,p_0)$ are the initial conditions)
B. The time evolution is governed by the Hamiltonian function $H(x,p, t)$.
According to the Hamilton-Jacobi theory, there exists Hamilton-Jacobi phase function $S(x_0, x_1, t_0, t_1)$ satisfying the Hamilton-Jacobi equation:
$$\displaystyle{\frac{\partial S}{\partial t}+H\left(x_1,\frac{\partial S}{\partial x_1}, t\right) = 0}$$
(where $(x_1,p_1)$ are the coordinates and momenta at time $t$)
The momenta can be derived from the Hamilton-Jacobi phase function:
$$\displaystyle{p_i = \frac{\partial S}{\partial x_i}}$$
The problem of expressing the state of the system in terms of the initial and final coordinates is rendered to a problem of transformation of probability distributions. We can define the state of the system in the initial and final coordinates as:
$$\displaystyle{F_t(x_0, x_1) = f\left(x_0,\frac{\partial S}{\partial x_1}(x_0, x_1, t) \right)}$$
The transformation Jacobian is given by:
$$ \displaystyle{dx_0 dp_0 = \frac{\partial^2 S}{\partial x_0\partial x_1}}dx_0 dx_1 $$
And the normalization condition:
$$\displaystyle{\int F_t(x_0, x_1) \frac{\partial^2 S}{\partial x_0\partial x_1}(x_0, x_1, t)dx_0dx_1 = 1}$$
In the general case, the joint distribution $F_t(x_0, x_1)$ will not be separable 
