  A calculation method for higher-order moments of physical quantities,
including magnetization and energy, based on the higher-order tensor
renormalization group is proposed. The physical observables are represented by
impurity tensors. A systematic summation scheme provides coarse-grained tensors
including multiple impurities. Our method is compared with the Monte Carlo
method on the two-dimensional Potts model. While the nature of the transition
of the $q$-state Potts model has been known for a long time owing to the
analytical arguments, a clear numerical confirmation has been difficult due to
extremely long correlation length in the weakly first-order transitions, e.g.,
for $q=5$. A jump of the Binder ratio precisely determines the transition
temperature. The finite-size scaling analysis provides critical exponents and
distinguishes the weakly first-order and the continuous transitions.
