Memristors have emerged as promising candidates for multilevel data storage, in-memory processing, and neural networks since their intrinsic programmability of resistance states under applied stimuli has been well revealed in memristor modeling. However, the programming uncertainty arising from the inherently stochastic nature of the device itself has been overlooked in previous modeling approaches. This omission hinders the incorporation of memristor stochasticity into time-domain circuit simulation. To address this issue, we propose a behavior model that incorporates real-time programming stochasticity. Our model stands out for several attributes: 1) programming stochasticity is included and exhibited in its resistance change over time; 2) its stochastic behavior is depicted by the summation of its deterministic behaviors and a noise signal; and 3) both deterministic behaviors and noise amplitudes depending on the pulse amplitude v and the memristor resistance R are determined by sufficient characterization data of our in-house TiO2 devices in a data-driven method. Consequently, our model is validated as highly matched to the characterized memristor device in terms of time-domain resistance evolution. Additionally, the modeling process can be adapted to different memristors with significant device variations. Furthermore, the model is transformed into the standard Verilog-A style for in-circuit simulation. To demonstrate its compatibility with system-level circuit simulation, a mixed-signal CMOS circuit is designed. This circuit explores the feasibility of storing multibit data within a single memristor, while considering its stochasticity.