Can I derive the Boltzmann distribution by an invariance argument? In statistical mechanics, the Boltzmann distribution gives the probability of a system being in state $i$ as
$$\displaystyle \frac{e^{- \beta E_i}}{\sum_i e^{-\beta E_i}}$$
where $E_i$ is the energy of state $i$.  I have generally seen this demonstrated, starting with some reasonable physical assumptions, via a heat bath argument (as exposited e.g. by Terence Tao) involving interactions between the system and a larger external system.  For me, an unsatisfying aspect of the heat bath argument is that it doesn't give me a strong reason to expect that a fundamental function like the exponential should appear at the end.
Here is what I think could be an argument which accomplishes that.  By inspection, the Boltzmann distribution only depends on the relative energies of the different states.  Under some mild assumptions this actually characterizes the Boltzmann distribution.  Let us suppose there is a non-negative function $f(E)$ such that WLOG $f(0) = 1$ and such that the probability of a system being in state $i$ is
$$\displaystyle \frac{f(E_i)}{\sum_i f(E_i)}.$$
Let us suppose that the system has two states.  Then the statement that the Boltzmann distribution only depends on the relative energies turns out to be equivalent to the functional equation $f(x + y) = f(x) f(y)$, which under any kind of continuity assumption whatsoever gives $f(x) = e^{ax}$ for some constant $a$. 
Question 1:  How can this argument be fleshed out?  In particular, what physical principle would suggest that the Boltzmann distribution only depends on the relative energies of the states?  (I seem to recall from my high-school physics lessons that energies are only well-defined up to an additive constant, but I would really appreciate some clarification on this issue.)
Question 2:  How does this argument relate to the heat bath argument or the combinatorial argument given, for example, at Wikipedia?  
(Motivation:  some important functions in mathematics, like the Jones polynomial and various zeta functions, can be interpreted as partition functions of certain statistical-mechanical systems, and I am trying to sharpen my physical intuition about these constructions.)
 A: Like Andreas, I find a maximum entropy argument to be intellectually appealing.  However, 
he says the solution can be found by Lagrange multipliers and I don't know the justification for using Lagrange multipliers. That is, in the space of all probability distributions on the particles, how do you know the maximum entropy solution is really accessible to variational methods?
For a derivation not using Lagrange multipliers, see the bottom of page 9 through page 11 
at http://www.math.uconn.edu/~kconrad/blurbs/analysis/entropypost.pdf.
A: For me the clearest derivation of the Boltzmann distribution is by maximizing the entropy $\sum n_i \ln(n_i)$ unter the constraint of constant total energy $\sum n_i E_i = \text{const.}$ and constant total particle number $\sum n_i = \text{const.}$. The Lagrange multiplicator for the first constraint gives $\beta$. You can immediately see that a shift of the energies does not change the distribution. 
A: I sketched this here and in section IV of my ancient preprint Effective statisical physics of Anosov systems.
A: This answer is just an expanding version of Kconrad answer's. I am posting it here because this argument support the variational method for finite state space and also touch in the observation made by Kconrad about a technical issue about boundary values of the variational approach. 

Proposition: Suppose $\Omega$ non empty finite set and let $\mathcal M$ denote the  set of the probability measures on $\Omega$ then
        $$
  \sup_{\mu\in\mathcal M} \left[ h(\mu)-\int_{\Omega} U d\mu \right]=\log Z
  $$
    moreover, the supremum is attained for the measure $\mu$ given by 
        $$
  \mu(\{\omega\})=\frac{1}{Z}e^{-U(\omega)}.
  $$
Proof: 
    Let be $n$ the cardinality of $\Omega$. Define the function $f:\mathbb R_+^n\to\mathbb R$ by 
        $$
  f(x_1,\ldots,x_n)=-\sum_{i=1}^n \Big[x_i\log x_i +K_ix_i\Big],
  $$
    where $K_i\in\mathbb R$ for all $i\in\{1,\ldots,n\}$. Consider the function $g:\mathbb R_+^n\to\mathbb R$ given by 
        $$
  g(x_1,\ldots,x_n)=\sum_{i=1}^n x_i.
  $$
    We fix an enumeration for $\Omega$ and let be $K_i=U(\omega_i)$. So the following optimization problem 
        $$
  \sup_{\mu\in\mathcal M} \left[ h(\mu)-\int_{\Omega} U d\mu \right]
  $$ 
    can be solved by finding a maximum for $f$ restricted to $g^{-1}(1)$. Note that for any critical point $(x_1,\ldots,x_n)$ of $f$ in $(0,\infty)^n\cap g^{-1}(1)$, it follows from the Lagrange Multipliers Theorem's that 
        $$
  \nabla f(x_1,\ldots,x_n)=\lambda \nabla g(x_1,\ldots,x_n)
  $$
    for some $\lambda\in\mathbb R$, i.e.,
        $$
  -(\log x_i +1+K_i)=\lambda, \ \ \ \text{for all}\ i=1,\ldots,n.
  $$
    So for any pairs of index $i,j\in\{1,\ldots,n\}$, we have 
        $$
  \log x_i +K_i=\log x_j+K_j
  $$
    taking the exponentials it follows that
        $$
      x_ie^{K_i}=x_je^{K_j}.
  $$
    Using that $\sum_{i=1}^nx_i=1$ and the above identities, we have 
$$x_ie^{-K_i}=\left[1-\sum_{j\in \{1,\ldots,n\}\backslash\{i\}}x_j\right]e^{-K_i}$$
$$=e^{-K_i}-\sum_{j\in \{1,\ldots,n\}\backslash\{i\}}x_je^{-K_i}$$
    So 
        $$
  x_ie^{-K_i}=e^{-K_i}-x_i\sum_{j\in \{1,\ldots,n\}\backslash\{i\}}e^{-K_j}.
  $$
    Explicting $x_i$, we show that all critical points of  $f$ in $(0,\infty)^n\cap g^{-1}(1)$ are given by (here there is just one)
        $$
  x_i=\frac{e^{-K_i}}{\sum_{j=1}^ne^{-K_j}}.
  $$
    The image of $f$ at this point is given by 
    $$
 -\sum_{i=1}^n \left[\left(\frac{e^{-K_i}}{\sum_{j=1}^ne^{-K_j}}\right)\log 
 \left(\frac{e^{-K_i}}{\sum_{j=1}^ne^{-K_j}}\right) 
 +K_i\left(\frac{e^{-K_i}}{\sum_{j=1}^ne^{-K_j}}\right)\right]
 =
 \log\left(\sum_{j=1}^ne^{-K_j}\right)
 $$
    to see that $(x_1,\ldots,x_n)$ is local maximum we can compute the Hessian and check that it is negative definite at this point.

To show that the point is global maximum point, we can compare the image of 
    $f$ at this point, with the value of  $f$ in any point of the set 
    $\partial (0,\infty)^n\cap g^{-1}(1)$. 
    The restriction of $f$ to this set is given by 
        $$
  f(x_1,\ldots,x_n)=-\sum_{i\in\{1,\ldots,n\}\backslash I}\Big[x_i\log x_i +K_ix_i\Big]
  $$
    Where $I\subset \{1,\ldots,n\}$ is a index set such that $|I|\geq 1$ e $x_i=0$ para todo
    $i\in I$ .
    We define $f_I$ which is a function of $n-|I|$ variables. It is maximum point can be determined in the same way and therefore we have that the max of $f_I$ is
$$
\log\left(\sum_{j\in\{1,\ldots,n\}\backslash I}e^{-K_j}\right)
$$
which is less than 
$$
\log\left(\sum_{j=1}^ne^{-K_j}\right).
$$ 
Repeating this argument at most $n$ times we conclude that maximum of $f$ restricted to $g^{-1}(1)\cap \mathbb R^n_+$, is not attained in the boundary.
A: I'd like to chime in here, as someone with a physics background.
I absolutely love the derivation given by Landau in volume 5 on statistical physics, chapter 1.  The basic idea is that since the log of the probability distribution function (i.e. the entropy) is an additive constant of the motion, it can be expressed as a linear combination of the 7 fundamental additive constants of the motion, namely the three components of momentum, the three components of angular momentum, and the energy.  But since the momentum/angular momentum components can be reduced to zero with an appropriate frame of reference, the log of the distribution function depends only on some multiple of the energy, which turns out to be 1/T.  We obtain the partition function naturally by normalizing the probability distribution.
I think this answers your question 1 from a physics point of view.
EDIT:
in view of the comments below, I should point out the the probability distribution I am referring to gives the probability of finding a system of N particles which obey the laws of classical mechanics in the state for which the nth particle is at position rn and moving with a velocity vn
A: This is an old question, but I would like to contribute a small note about energies being defined up to a constant, from a different point of view.
Force is usually seen as a vector, but you can see it as a 1-form, which integrates along a curve to produce the work done. A conservative force is one for which work is independent of path, i.e. it integrates to zero along any closed curve. By Stoke’s theorem, its exterior derivative vanishes  $dF=0$, i.e. it is a closed form. By the Poincaré lemma, it is exact, $F=dU$, where $U$ is the energy. Shift invariance, a kind of Gauge invariance, is now obvious.
