  In this work, we develop a stochastic matrix product state (stoMPS) approach
that combines the MPS technique and Monte Carlo samplings and can be applied to
simulate quantum lattice models down to low temperature. In particular, we
exploit a procedure to unbiasedly sample the local tensors in the matrix
product states, which has one physical index of dimension $d$ and two geometric
indices of dimension $D$, and find the results can be continuously improved by
enlarging $D$. We benchmark the methods on small system sizes and then compare
the results to those obtained with minimally entangled typical thermal states,
finding that stoMPS has overall better performance with finite $D$. We further
exploit the MPS sampling to simulate long spin chains, as well as the
triangular and square lattices with cylinder circumference $W$ up to 4. Our
results showcase the accuracy and effectiveness of stochastic tensor networks
in finite-temperature simulations.
