Can energy constraint $\sum n_r\varepsilon_r=E$ for nearly independent particle systems be replaced by $1/E$ or any nonlinear function $f(E)$? Background Information
For a nearly independent particle system, $N$ denotes the number of particles, $E$ the total energy. Suppose there are $n_r$ particles at a specific energy level $\varepsilon_r$, then according to particle and energy conservation,
$$
\sum_r n_r=N,\quad \sum_r n_r\varepsilon_r=E.
$$
These two equations are constraints for the system. Any sequence of $n$'s forms a set $\{n_r\}$, called a distribution. There can be many microstates corresponding with a distribution, the number denoted by $W(\{n_r\})$. For identical localized particles,
$$
W(\{n_r\})=N!\prod_r\frac{g_r^{n_r}}{n_r!},
$$
where $g_r$ denotes the degeneracy of energy level $\varepsilon_r$.
To determine the most probable distribution, Lagrange multiplier method is applied and gives
$$
\delta W(\{n_r\})-\alpha \delta\sum_r n_r-\beta\delta \sum_rn_r\varepsilon_r=0.
$$
With several approximations ($n_r\gg1$ and Stirling's formula), we end up with
$$
n_r\propto\exp(-\beta \varepsilon_r),
$$
or
$$
P(\varepsilon_r)\propto\exp(-\beta \varepsilon_r),
$$
where $P(\varepsilon_r)$ denotes the probability of a particle at energy $\varepsilon_r$.

Question
Now my question is, can the second constraint, $\sum n_r \varepsilon_r=E$, be replaced by other constraints such as $\sum n_r f(\varepsilon_r)=f(E)$ where $f(E)$ is a non-linear function?
The new constraint doesn't require additivity of $f(E)$, since the second constraint can be regarded as the average of some physical quantity. In the original form, the quantity is total energy $E=N\varepsilon$, where $\varepsilon$ is the average, so
$$
\varepsilon = \sum_r\frac{n_r}{N}\varepsilon_r,
$$
with subsequent derivation unchanged.
In the new form, the quantity is chosen to be $f(E)$, whatever it is, for example $f(E)=1/E$. And the average is $f(E)/N$. So
$$
\frac{f(E)}{N}=\sum_r\frac{n_r}{N}f(\varepsilon_r) \Rightarrow f(E)=\sum_r n_r f(\varepsilon_r).
$$
Of course for a given distribution, the average $1/E$ is not the reciprocal of the average $E$.
But once $E$ is fixed, so is $f(E)$. No matter how the distribution $\{n_r\}$ is changed,
$f(E)$ conserves.
The same derivation remains legal and gives
$$
P(f(\varepsilon_r))\propto \exp(-\beta' f(\varepsilon_r)).
$$
These two probability is certainly not the same, since
$$
P(\varepsilon_r)\mathrm d\varepsilon_r=
P(f(\varepsilon_r))\mathrm df(\varepsilon_r)
\Rightarrow
P(\varepsilon_r)=
P(f(\varepsilon_r))f'(\varepsilon_r).
$$
What's wrong? Thanks for answering.

In Pathria's statistical Mechanics, section 3.2, he uses the same trick to derive the distribution of canonical ensembles. In an ensemble, there's no "total energy" but an average one, and it seems uncessary (although natural) to prefer $E$ rather than $f(E)$, which draws me to this question.
 A: It might be possible to proceed like this. I just see this procedure as a change of variables. Look at it the other way around: given that energies are distributed according to
$$P(\varepsilon)=Z^{-1}g(\varepsilon)e^{-\beta \varepsilon},$$
You can compute the probability distribution of any other variable $f=f(\varepsilon)$ as
$$P(f)=P(\varepsilon(f))\frac{d\varepsilon}{df}.$$
In any case, if you propose to fix $f(\epsilon)$, then you would need to look for the "thermodynamic meaning" of $\beta'$. Consider that $\beta'$ might no longer be the inverse of the temperature.
A: My point now is, the energy constrained cannot be replaced by another one.
Because the sytem is isolated, its energy
$$
\langle E\rangle=\sum n_r\varepsilon_r
$$
is restricted to be constant.
This means, we allow the system to transition between microstates of equal energy.
If we replace the constraint with a non-linear function $f(E)$,
what we restrict is
$$
\langle f(E)\rangle=\sum n_rf(\varepsilon_r).
$$
And the system is allowed to transition between microstates that conserve $f(E)$.
Mathematically, these two constraints are not equivalent, as stated in the original post.
When the system transition between microstates of equal $f(E)$, the total (or average) energy of the system is not conserved, since $f(E)$ is non-linear.
Physically, it's not what we want for an isolated system.

The problem with nearly independent particle systems can be explained this way. In Pathria's book, though, there's no need to isolate the whole ensemble, and the problem still exists.
