Simulating thermal-equilibrium properties at finite temperature plays a crucial role in studying the nature of quantum many-body systems. In particular, implementing a finite-temperature simulation on a quantum computer is expected to overcome the difficulty in simulating large-sized systems, for which the quantum Monte-Carlo technique on a classical computer suffers from the sign problem in general. While several methods suitable for fault-tolerant quantum computing (FTQC) devices are expected to be useful in studying large-scale quantum many-body systems, those proposed so far involve a large number of ancilla qubits and a deep quantum circuit with many basic gates, making them unsuitable for the early-FTQC era, i.e., the early stage of the FTQC era, at which only a limited number of qubits and quantum gates are available. In this paper, we propose a method suitable for quantum devices in this early stage to calculate the thermal-equilibrium expectation value of an observable at finite temperature. Our proposal, named the Markov-chain Monte Carlo with sampled-pairs of unitaries (MCMC-SPU) algorithm, is based on sampling simple quantum circuits and generating the corresponding statistical ensembles, and overcomes the difficulties in the resource requirements and the decay in probability associated with postselection of measurement outcomes on ancilla qubits. We demonstrate the validity of our proposal by performing a numerical simulation of the MCMC-SPU algorithm on the one-dimensional transverse-field Ising model as an illustrative example.