Metal halide perovskites (MHPs) have soft lattices with strong anharmonicity and will undergo entropy-driven solid−solid phase transitions upon heating. Here, we investigate the polymorph stabilities and phase transitions in one of the lead-free MHPs, CsSnI 3 , by several molecular simulation techniques. Three different phase transitions (γ ↔ β, β ↔ α, and yellow → black) in CsSnI 3 have been successfully reproduced by molecular dynamics (MD) simulations with a newly developed empirical force field. The heating and annealing MD simulations and free-energy calculations with the non-equilibrium thermodynamic integration (NETI) method predict the transition temperatures of 275, 385, and 280 K for the γ ↔ β, β ↔ α, and yellow → black transitions, respectively. Lattice dynamics (LD) simulations within the harmonic approximation fail to predict the correct phase stability in CsSnI 3 at high temperatures. The quasiharmonic approximation (QHA) calculations that include the volume dependence of the phonon frequencies and lattice energies correctly predict all phase transitions in CsSnI 3 . However, the transition temperatures of the γ ↔ β and β ↔ α transitions predicted by the QHA calculations significantly deviate from those by MD simulations. By comparing the Gibbs free energies calculated by the LD simulations within the QHA and MD-based NETI method, we find the differences of 3−30 meV for different polymorphs. Although calculations based on the harmonic model can provide valuable information, the anharmonic terms need to be included for accurate predictions of transition temperatures of phase transitions in CsSnI 3 and other MHPs.