{'fraction_non_alphanumeric': 0.03281267248040622, 'fraction_numerical': 0.010652389888508666, 'mean_word_length': 4.414973102211596, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Extended Kramers-Moyal analysis applied to optical trapping

The Kramers-Moyal analysis is a well established approach to analyze stochastic time series from complex systems. If the sampling interval of a measured time series is too low, systematic errors occur in the analysis results. These errors are labeled as finite time effects in the literature. In the present article, we present some new insights about these effects and discuss the limitations of a previously published method to estimate Kramers-Moyal coefficients at the presence of finite time effects. To increase the reliability of this method and to avoid misinterpretations, we extend it by the computation of error estimates for estimated parameters using a Monte Carlo error propagation technique. Finally, the extended method is applied to a data set of an optical trapping experiment yielding estimations of the forces acting on a Brownian particle trapped by optical tweezers. We find an increased Markov-Einstein time scale of the order of the relaxation time of the process which can be traced back to memory effects caused by the interaction of the particle and the fluid. Above the Markov-Einstein time scale, the process can be very well described by the classical overdamped Markov model for Brownian motion.


I. INTRODUCTION
Many real world stochastic processes q(t) can be modeled by stochastic differential equations of Langevin type [1][2][3],q = D (1) (q) + 2D (2) where D (1) (q) is called the drift coefficient that describes deterministic influences on the dynamics and D (2) (q) is the diffusion coefficient that describes the state dependent amplitude of a fast fluctuating force Γ (t) with zero mean, adopting Itô's interpretation. A frequently applied idealization is to assume that Γ (t) is δ correlated in time and Gaussian distributed. In this case the process fulfills the Markov property and can be completely characterized by the corresponding Fokker-Planck equation (FPE) with the Fokker-Planck operator that describes the temporal evolution of the probability density function (PDF) f q (x, t) = δ(x − q(t)) of the process q(t). The Fokker-Planck operatorL(x) depends on the drift and diffusion coefficients that determine the corresponding Langevin equation (1). Drift and diffusion coefficients are also referred to as the first and second Kramers-Moyal (KM) coefficients, * c.honisch@uni-muenster.de respectively. Given the transition probability densities p q (x ′ , t + τ |x, t) of a general Markov process q(t), the nth KM coefficient can be defined as where M (n) τ (x) is called the nth conditional moment given by τ (x) = (q(t + τ ) − q(t)) n | q(t)=x = (x ′ − x) n p q (x ′ , t + τ |x, t)dx ′ .
Since the transition PDFs can be estimated from measured time series data of a process of interest, it is in principle possible to set up a model in terms of a Fokker-Planck or Langevin equation by data analysis. Since Friedrich and Peinke applied this nowadays called KM analysis to the investigation of the turbulent cascade [4,5] in 1997, it has developed to a rapidly growing field of research with many applications in the natural sciences and beyond (see [6] and references therein). One major problem connected to the KM analysis is the limit of the time increment τ to zero that has to be performed in the determination of the KM coefficients. Without an appropriate limiting procedure, estimated KM coefficients can deviate significantly from the true coefficients if the minimal available τ , i. e. the sampling interval of the measured time series, is too large. These errors are referred to as finite time effects and were the subject of many publications during the last years (e. g. [7][8][9][10][11]). Without an adequate understanding of finite time effects, KM analysis involves the risk of dramatic misinterpretations of the achieved results.
In a recent publication [12], we have developed a method that allows for a correct KM analysis when the sampling interval is large. In the present article we demonstrate that also this method fails if the sampling interval is so large that the information about the true KM coefficients is (almost) lost. In order to lower the risk of misinterpretations in those cases, we extend the method by the computation of error estimates for the obtained model parameters. The problem of error estimates for model parameters in the KM analysis context is also addressed in a very recent publication by Kleinhans [13] in a Bayesian framework. In contrast to our approach, this method is only designed for data sets with a sufficiently large sampling rate.
We present an application of our method to an optical experiment that yields real-world stochastic data via trajectories of a Brownian particle trapped by optical tweezers. We apply our method to these trajectories in order to yield estimations for the forces induced by the optical tweezers and to analyze the spatial distribution of the temperature.
The outline of the article is as follows. In Sec. II we review the problem of finite time effects and discuss an example in order to provide an intuitive understanding. Afterwards we review the method of Ref. [12] and present important extensions. Finally, Sec. IV describes the application to our experiment. The last section is dedicated to some concluding remarks.

II. FINITE TIME EFFECTS
At first we introduce the finite time KM coefficients where M (n) τ (x) are the conditional moments defined in Eq. (5). Apparently, lim τ →0 D (n) τ (x) = D (n) (x). In order to understand finite time effects and to possibly correct for them, it is important to know, how the finite time coefficients change with the time increment τ given the true coefficients.
There are several possibilities to compute finite time coefficients. One possibility is by solving the FPE. The transition PDF p q (x, t 0 + τ |x 0 , t 0 ) that occurs in Eq. (5) is the solution to the FPE (2) at time t = t 0 + τ with the initial condition Unfortunately this is in most cases impossible to do analytically. Also for a numerical solution of the FPE, the initial condition in form of a Dirac δ distribution will cause problems. Another possibility is to numerically integrate the corresponding Langevin equation and to estimate the finite time coefficients from the generated time series. Based on this approach, Kleinhans et al. developed a maximum likelihood approach to estimate KM coefficients for processes with finite sampling rates [14,15].
A third possibility is the adjoint operator approach. In Ref. [8], Friedrich et al. mention that the conditional moments can be expressed as a Taylor series in τ which reads This series expansion contains the adjoint Fokker-Planck With the adjoint operator method it is not only possible to obtain a series expansion. As Lade showed in Ref. [9], the conditional moments can also be computed by solving the partial differential equation with initial condition Then the conditional moments are given by Eq. (10), also called the adjoint FPE (AFPE), can easily be solved analytically as long as the drift is linear. Otherwise numerical solutions can be obtained with standard finite difference schemes.
For an Ornstein-Uhlenbeck (OU) process, the analytic expressions for the finite time coefficients read We have plotted the finite time coefficients together with the conditional moments in Fig. 1 (drift) and Fig.  2 (diffusion) for γ = α = 1. By taking a logarithmic scaling for the τ axis, one can clearly identify two limiting cases separated by the relaxation time, which is τ R = 1 in this case. For τ ≪ τ R , the finite time KM coefficients have converged to the true KM coefficients, whereas the conditional moments tend to zero. For τ ≫ τ R , the conditional moments become stationary with respect to τ , whereas the corresponding KM coefficients vanish. The latter case is the limit of statistical independence that has been discussed by Anteneodo and Queirós in Ref. [10]. For τ ≫ τ R , the realization x of a process q at time t becomes statistically independent from the realization x ′ at time t + τ . Therefore, For a stationary process, f q (x ′ , t+τ ) is just the stationary PDF f q (x). As a consequence, the conditional moments yield and the finite time KM coefficients read That means that every Langevin process, independent of the true KM coefficients, appears as a process with linear drift and quadratic diffusion if the sampling interval of the data set is large compared to the relaxation time of the process. Therefore, as a preinvestigation before a KM analysis, it is necessary to investigate the autocorrelation function of a process in order to estimate the relaxation time, and to estimate the conditional moments for various available τ and present it as a plot like Figures 1 and 2. Furthermore, Ref. [11] provides some checks in order to validate models with linear drift and quadratic diffusion.
We can conclude that in the regime τ ≪ τ R , finite time effects are small and low order corrections should suffice or can even be neglected. In the regime τ ≫ τ R , all information about the true dynamics is lost and a KM analysis is no longer possible. In between those two regimes, if τ ≈ τ R , a KM analysis is possible with use of our method presented in Ref. [12], as is demonstrated in several examples. Of course, at some point depending on the available amount of data, the results obtained by this method will also fail when τ approaches the limiting case of statistical independence. The point at which it fails will depend on the size of the available data set. The lower statistical uncertainties are in the estimation of the conditional moments, the more accurate one can extrapolate to τ = 0. To quantify this possible accuracy, we extend our method by the computation of errors in the estimated parameters with which drift and diffusion coefficients are parametrized. This is done by a Monte Carlo error propagation technique that is often used in nonlinear optimization problems [16].
In the next section, we first summarize the general method and present some minor improvements. Afterwards the Monte Carlo error propagation technique will be introduced in some more detail. We also address the limitations of the method due to finite amount of data and large sampling intervals.

III. KM ANALYSIS AT LOW SAMPLING RATES
A. Review of the method and minor improvements The optimization procedure introduced in [12] is based on an inversion of Lade's method [9] to compute finite time coefficients with help of the AFPE, which was already mentioned in the foregoing section. The basic steps that have to be performed are the following.
As a first step, one has to obtain estimations of the conditional momentsM (1,2) τ (x) for several values of x and τ from the time series data. For this we recommend a kernel based regression with use of the Nadaraya-Watson estimator. The bandwidth was in [12] selected according to Silverman's rule of thumb [17] as whereσ is the standard deviation of the time series data. We found that this bandwidth is in many cases a bit too large leading to underestimated slopes of straight lines or curvatures of parabolas. More sophisticated data driven bandwidth selectors such as cross validation techniques or the method described in [18] are computationally expensive and did not yield robust improvements in our tested cases. We obtained the best results by simply reducing the factor 1.06 in (18) to 0.8. The next step is to find a suitable parametrization for the drift and diffusion coefficients. In most cases a polynomial ansatz will be appropriate. Also a representation in form of spline curves has been successfully tested in [12], where the interpolation points serve as optimization parameters. In either case, the coefficients are represented as D (1,2) (x, σ) with a set of optimization parameters σ. If an analytic expression has been chosen for the parametrization, one can try to obtain an analytic expression for the corresponding conditional moments M (1,2) τ (x, σ) via Lade's method. Otherwise, they have to be calculated by numerical integration of the AFPE. The optimal set of parameters must then be obtained by minimizing a distance measure between the estimated conditional momentsM (1,2) τ (x) and the computed ones. In [12] we have chosen the weighted least squares distance whereσ (1,2) ij are the statistical errors in the estimation of the conditional moments. In [12] we used as an error estimate, where T is the number of data points of the time series data, K(•) is the kernel function, and h is the selected bandwidth. The error estimate (20) does not take into account the influence of the bandwidth. According to [17], Eqs. (20) have to be corrected by a factor of where ||K|| 2 2 is the L 2 norm of the kernel, thus For the Epanechnikov kernel which we use, ||K|| 2 ≈ 0.518. However, regarding the optimization, the correction factor (21) has no influence, since the error estimate has only the purpose of a relative weighting of the estimated conditional moments for different values of x and τ . It only becomes relevant in connection with the Monte Carlo error propagation described in Sec. III B. The minimum of (19) is in [12] determined by a trust region algorithm, but other optimization methods should work as well.

B. Monte Carlo error propagation
The basic idea of the Monte Carlo error propagation (MCEP) approach is as follows: After estimating the model parameters of a specific model from a noisy data set (backward problem), the model is used to produce an unnoisy data set (forward problem). Then one generates an ensemble of pseudo data sets by perturbing the unnoisy data set with different realizations of noise. For each data set again, the model parameters are estimated. Eventually, one can compute the standard deviation for each model parameter from this ensemble, which can be used as an uncertainty measure.
In the present case, we do not regard the time series data as our noisy data set, but the estimated conditional momentsM (1,2) τ (x). In order to obtain reliable uncertainty estimates with the MCEP approach, it is crucial that the artificial perturbations of the unnoisy data set are of the same amplitude as the errors of the original data set. Therefore one needs a reliable absolute error estimate for the conditional moments and the correction factor (21) becomes important.
To generate the ensemble of pseudo data sets, one has to compute the conditional moments that correspond to the estimated model parameters by solving the corresponding AFPE. To each M (1,2) τi (x j ) we add an independent Gaussian distributed random number with standard deviation of expression (20) multiplied by the factor (21) for each member of the data set ensemble. Then the model parameters are estimated again for each data set and the corresponding standard deviation is calculated for each model parameter.
There is one technical problem connected to the MCEP approach. It can happen that the optimization routine does not find the absolute minimum of the least square potential for some pseudo data sets of the ensemble. Although the trust region algorithm is a very robust one, our observation is that it sometimes gets stuck in local minima close to the initial condition. If the initial condition is selected the same for all pseudo data sets, it can therefore happen that the computed standard deviation is significantly underestimated. To overcome this problem, we choose random initial conditions that lie in a region around the estimated parameters from the original data set. This in return bears the danger that the optimization routine will not find the way back to the true minimum if the initial condition lies too far away. This can lead to overestimated standard deviations. Therefore, only optimization results are accepted where the final residual is below a specific threshold. Otherwise the optimization is repeated with another initial condition. The threshold is selected three times the residual of the result from the original data set. Nevertheless it can happen that some outliers lead to a significantly overestimated standard deviation, especially in cases of very large sampling intervals. The easiest way to overcome this problem is to detect these outliers by eye and remove them from the ensemble. To test the MCEP approach, we generate an ensemble of ten synthetic time series of an OU process with D is stored so that the sampling interval of time series is τ = 1. Each time series consists of 10 6 data points. For the optimization we choose the parametrization with the optimization parameters γ, α and β. For each of the ten data sets, we estimate these parameters as well as the corresponding standard deviations with the MCEP approach. The results are depicted in Fig. 3. By eye, the estimated error bars seem reasonable. To make a quantitative statement, we compute the standard deviations of the estimated parameters from the ensemble and compare them to the mean standard deviations of the MCEP approach. The numbers are listed in Table I. The errors of the parameters γ and α fit very well while the error of β is a bit overestimated. However, one should note that the errors obtained by the MCEP approach only cover uncertainties caused by statistical fluctuations due to the finite amount of data. The influence of other error sources such as measurement noise, deviations from the Markov property, non-stationarity of the time series or an inappropriate parametrization of drift and diffusion coefficients, just to name a few, are not considered. Therefore, the MCEP errors should be regarded as a lower bound for the true errors. Nevertheless they can decrease the danger of an overestimation of the significance of KM analysis results.

C. Limitations of the approach
Now we use the MCEP approach to demonstrate the limitations of the KM analysis caused by finite time effects and limited amount of data. To this end, we compute the error estimates for the parameters γ, α and β for synthetic time series with different sampling intervals. The synthetic time series are generated in the same manner as described in section III B. All data sets consist of 10 6 data points. The result is shown in Fig. 4. One can see that the estimation of the diffusion parameters α and β breaks down around 2.5 relaxation times (τ R = 1), while the drift parameter γ can still be estimated up to about 5 relaxation times. Above these values it is difficult to obtain robust error estimates with the MCEP approach because the afore-mentioned outliers become very frequent. However, using the MCEP approach for sampling intervals above these values, one will clearly notice that the KM analysis is no longer feasible.
Looking at Fig. 4, one further notices that the errors of α and β are strongly correlated, which is not surprising. That the absolute sizes of the errors are almost equal is only the case for the specific selection of parameters γ = α = 1 and β = 0.

A. Experimental setup
Trapping of small dielectric objects in the nano-and microscale can be realized by focusing a single laser beam so that it creates a gradient force near the focal region that is able to hold these particles, and even levitate them against gravity [19]. Our trapping experiment is based on such a single beam gradient laser tweezers system (see FIG. 5) as in [20]. with a Gaussian beam distribution is coupled into an inverted microscope (Ti Eclipse, Nikon) for trapping with an oil-immersion objective (CFI Apo TIRF 100x) with a numerical aperture of NA = 1.49. This objective is also used for observation purposes and images the observation plane onto a high-speed CMOS sensor (MV2-D1280-640-CL-8, Photonfocus) with a frame rate of 488 Hz at full resolution. By reducing the size of the captured images to a region of interest of 80 x 64 pixels, the framerate could be increased to 3966 Hz, which corresponds to a sampling interval τ s ≈ 0.25 ms. For the measurements uncoated polystyrene beads (Kisker) with a diameter of (1.002 ± 0.043) µm are used. These microparticles are suspended in destilled water and trapped 5-10 µm above the cover glass. Due to technical reasons up to 150 videos consisting of 10 4 frames are captured over 15 minutes and analyzed afterwards. The particle's position is determined by a Matlab algorithm implementing a center of mass detection [21,22].
The autocorrelation function of a trapped particle's x(y)-position decays with a relaxation time τ x(y) = 6πηr/k x(y) where η is the viscosity of the surrounding fluid, r the radius of the particle and k x(y) the corresponding trap stiffness [23]. For a particle with a radius of r = 0.5 µm typical relaxation times are around 10 ms µm pN −1 /k x . Therefore the laser power is reduced to about 10 mW in the trapping plane which results in a trapping stiffness of k x = 21 pN/µm and thus a relaxation time of about τ x = 0.45 ms (compare Fig. 9) in order to have sufficient temporal resolution to observe the particle's dynamics.

B. Modeling
The equation of motion for a Brownian particle in an optical trap reads for one spatial dimension [25] where m p is the mass of the particle, F fr is the friction force, F op is the force induced by the optical trap and F th denotes the thermal fluctuations. The generell form of the friction force reads According to the fluctuation dissipation theorem [24], the kernel γ(t) is connected to the correlation function of the thermal fluctuations Assuming a laminar velocity profile around a spherical particle and no-slip boundary conditions, the kernel is given by where η is the dynamic viscosity of the fluid and r the radius of the particle. This yields the well-known Stokes' law with λ = 6πηr. On time scales τ ≫ τ I = m p /λ (τ I ∼ 0.05 µs in our experimental setup), inertia can be neglected. With a linear optical force F op = −kx, this yieldsẋ with Γ (t)Γ (t ′ ) = δ(t − t ′ ) and the diffusion coefficient which is known as the Einstein-Stokes equation. For a constant temperature equal to the room temperature of T = (294 ± 2) K, a particle diameter of 2r = (1.002 ± 0.043) µm, and a viscosity of η =(1.00±0.02) · 10 −3 Nsm −2 , the Einstein-Stokes equation predicts a diffusion constant of A more realistic treatment goes beyond Stokes' law and also considers the momentum that is transferred from the particle to the fluid. Ref. [25] gives a derivation of Eqs. (25), (26), and (27) for a macroscopic sphere in an incompressible fluid from linearized stochastic hydrodynamic equations. Thereby the Fourier transform of the memory kernel is computed aŝ Here, ρ f denotes the mass density of the fluid. This corresponds to a friction force where m f = (4π/3)ρ f r 3 is the mass of the displaced fluid. The correlation of the thermal fluctuations shows a negative algebraically decaying tail for |t − t ′ | −3/2 > 0. The memory term in the friction force (34) and the long correlations of the thermal fluctuations clearly obliterate the Markov property of the process. However, it is possible that a finite Markov-Einstein (ME) time scale τ ME exists, which means that the process can approximately be described by a Markov process for times larger than τ ME . If such a time scale is found, one can try to find a model in form of Eq. (30) with an effective drag coefficient λ and diffusion coefficient D (2) with our data analysis method.
To obtain a rough estimation of the ME time scale, we compare the mean squared displacement M SD(t) = (∆x(t)) 2 of the full hydrodynamic model, which was computed by Clercx [26], with the MSD of the Markov model, Eq. (30). Fig. 6 shows a plot of the MSDs for both models with the parameters according to our experiment. In the overdamped Markov model (red  In the full hydrodynamic model the smallest characteristic time scale is the inertia time scale τ I ∼ 10 −8 s. For times smaller than τ I the MSD grows quadratically with time, M SD H (t ≪ τ I ) = (k B T /m * )t 2 , where m * = m p + m f /2. This is the so called ballistic regime [27]. Above the ballistic regime it takes about four decades until the MSD approaches the diffusive regime and coincides with the overdamped Markov model. These four decades are influenced by hydrodynamic memory effects. The inset of Fig. 6 shows the relative deviation between the MSDs of the two models. According to this, the influence of the hydrodynamic memory effects is still present at times t ∼ τ k . For times larger than 0.4 ms, the relative deviation becomes smaller than one percent. Therefore, we expect a ME time scale of this order of magnitude. A direct test of our data for Markovianity is presented in Sec. IV D yielding a ME time scale of τ ME ≈ 0.5 ms which is in good agreement to the discussion above.
C. Data preparation Fig. 7 (c) shows the first twelve seconds of the x component of the original data set. One notices that the mean value of the time series fluctuates over time, i. e., the time series is not stationary. These fluctuations of the mean value lead to large values in the power spectral density (PSD) of the time series at low wave numbers k as one can see in Fig. 7 (a), which shows the PSD averaged over the 150 measured trajectories. To enhance the stationarity of the time series, the corresponding Fourier coefficients for each of the 150 trajectories are reduced such that they fit in the Lorentz-like form of the spectrum. The zero wave number coefficients are set to zero to center time series around zero. This kind of high pass filtering leads to the averaged power spectral density depicted in Fig. 7 (b). The corresponding time series after this preprocessing is shown in Fig. 7 (d). The same type of preprocessing is applied to the y component of the trajectories.

D. Preinvestigations
As a first preinvestigation, we take a look at the Markov property. A necessary condition for a process to be Markovian on a specific time scale τ is the validity of the Chapman-Kolmogorov equation (CKE) [1]: where p τ (x ′ |x) := p(x ′ , t + τ |x, t) denotes the transition probability density function of the process. We further To test our data for Markovianity on a time scale τ , we compare estimates for the joint PDFs f 2τ and f CK 2τ (cf. Ref. [6]). Fig. 8 shows the corresponding contour plots for τ = τ s and τ = 2τ s for the two components of the process. For τ = τ s (upper panels), one can see clear de-viations between the contour lines which indicate that the process is not Markovian on this time scale. For τ = 2τ s the deviations vanish. This leads to the conclusion that the process has a ME time scale τ ME ≈ 2τ s . Therefore we only include conditional moments with time increments τ ≥ 2τ s into the minimization of Eq. (19).
As a next step, we take a look at the autocorrelation functions (ACF) to get an impression of the typical time scales of the system and to decide whether the sampling interval τ s is sufficiently small for a reliable KM analysis. The ACFs are shown in Fig. 9. To evaluate the relaxation time as a typical time scale of the system, we fit an exponential, e −cτ , to the first points of ACF and take τ R ≈ c −1 as a rough estimate. According to this estimate, τ R ≈ 2τ s . This means that the relaxation time of the process is approximately equal to the ME time scale. Therefore, according to the discussion in Sec. III C, the KM analysis should be possible. As a last step, we check whether it is possible or not to regard the x and y components of the particle mo- tion as two independent processes. To this end, we measure the finite time drift vector field D (1) (x, y, τ ) = (D (1) x (x, y, τ ), D (1) y (x, y, τ )) ⊤ with for τ = τ s . We assume that the qualitative form of the drift vector field is not affected that much from finite time effects and deviations from the Markov property.
In Fig. 10 we show contour plots of D (1) x (x, y, τ ) (left) and D (1) y (x, y, τ ) (right). In a region with a radius of about 0.04 µm around the origin of the coordinate system, where most of the data points are, the contour lines of D (1) x (x, y, τ ) and D (1) y (x, y, τ ) are approximately linear and parallel to the y axis and x axis, respectively. Hence, Therefore, we treat the two components of the particle motion as independent processes.

E. Analysis results
To estimate the drift and diffusion coefficients for each of the two processes, we first make a parametric ansatz in form of an OU process, Eqs. (13), with optimization parameters γ and α. Our experience from synthetic time series data is that the best results are achieved if one includes the finite time coefficients of the n smallest time increments τ i = iτ s , such that τ n is between one and two relaxation times. Since we have a finite ME time scale of 2τ s in our case, we include the conditional moments with τ = iτ s , i = 2, 3, 4 into the least squares potential, Eq. (19).  Table II shows the results of the optimization as well as the estimated errors by the MCEP method for the x and y components of the process. A graphical representation is depicted in Figures 11 and 12, respectively. A possible spatial dependence of the temperature due to the heating of the laser cannot be resolved on the basis of the experimental data. If one includes a quadratic term in the diffusion ansatz as in Eq. (24b), the error estimate for the parameter β is of the same size as the estimated value. However, our model with linear drift and constant diffusion describes the process very well, as we will see in Sec. IV G.
The diffusion coefficients deviate by a factor of approximately 1.3 from the result that was expected from the Stokes-Einstein equation (s. Sec. IV B). To understand this deviation, we also measure the diffusion coefficients of different freely diffusing particles. The results are presented in the following section. The stiffness of the trap can nevertheless be calculated by If we assume a temperature of (294±2) K, we obtain In order to understand the deviations from the Einstein-Stokes equation, we also determine the diffusion coefficients for freely diffusing particles. In the same manner as described in Sec. IV A, we measure the positions of seven different particles that are not trapped by optical tweezers. For each particle, one time series with approximately 10 4 time steps at a sampling frequency of 3873 Hz is measured. From this data we determine the mean squared displacement (MSD) for which the relation holds. Fig. 13 shows the obtained MSDs for the x components of the seven particles. The diffusion coefficients can which is in good agreement to the expected result according to the Einstein-Stokes equation, Eq. (32). The reason for the higher standard deviation of 0.06 (µm) 2 s −1 might be that the fluctuations among the particle radii are larger than indicated by the manufacturer. However, the low diffusion coefficient found in the optical trapping experiment is inside the range of fluctuations of diffusion coefficients among different particles.

G. Comparison between model and data
In this section we compare our estimated models to the experimental data. At first we compare the single point PDFs. The PDFs of the experimental trajectories are estimated via a standard kernel method using the Epanechnikov kernel (23) and a bandwidth according to Eq. (18). The PDFs according to the OU model are given by ] t [s] particle 1 particle 2 particle 3 particle 4 particle 5 particle 6 particle 7 Einstein-Stokes time series for the two processes. In both cases the PDFs agree very well. The second quantity we compare is the ACF, Eqs. (39). Fig. 15 shows the ACF of the experimental time series for the x and y components of the particle motion, together with ACFs of our model that are given by As one can also see in Fig. 9, the experimental ACFs are not exactly exponentially shaped. The reason for this is discussed in Sec. V. However, one can see that the relaxation times of the models and the experimental data sets coincide very well. As a next step, we compare the first and second finite time conditional moments. They are depicted in Figs  are qualitatively equal. The conditional moments of our model are given by Apparently, our model fits very well to the data. Significant deviations do only occur for the smallest τ = τ s , which is below the ME time scale and was not included in the optimization. Here the slope (with respect to x) of the first moment and the second moment (for all x) of the data are smaller than in the model. An inclusion of time increments below the ME time scale would therefore lead to underestimated drift and diffusion coefficients. Instead of the fourth conditional moments, we compare the fourth finite time KM coefficients which reads for our model D (4) τ (x) = 1 24τ (1 − z τ ) 4 x 4 Here we have used the abbreviation z τ = e −γτ . Fig.  18 shows the estimated fourth KM coefficient together with Eq. (49). As in the two previous plots, significant deviations are only visible for the smallest τ = τ s , where the coefficient estimated from data is smaller than in the model. However, one can see that D (4) τ (x) vanishes as τ approaches zero. This is the requirement for the Pawula theorem that guarantees that also the third and all higher KM coefficients vanish [1].

V. CONCLUSION
In the present article we have demonstrated that finite time effects can lead to significant quantitative and also qualitative errors in the KM analysis. A previously published method [12] allows in principle to correct for errors caused by finite time effects. But even the application of this method bears the danger of misinterpretation of the achieved results if the sampling interval of an experimental time series data set approaches the limit of statistical independence discussed in [10]. In order to avoid these misinterpretations, we have extended this method by the calculation of error estimates for the determined model parameters. These error estimates allow for a more honest assessment of the validity of the obtained results. Since the estimated errors only take into account errors caused by finite time effects and the finite amount of available data, they should be regarded as lower bounds for the true model errors, which can also be influenced by other effects.
As long as analytic solutions of the AFPE for the selected parametrization of the KM coefficients are known, the computational effort is very low. The applications presented in this articles have a computation time of less than five seconds on a usual desktop computer. If the AFPE has to be solved numerically, the computational effort increases dramatically. The MCEP method has not been tested for those cases.
To apply our method to real-world stochastic data, we have also conducted an experiment where trajectories of Brownian particles trapped by an optical tweezers system were measured. We find a ME time scale which is approximately equal to the relaxation time of the process. This can be explained by hydrodynamic memory effects that are still present on this time scale. The large ME time has the consequence that even if the trajectories of the particle were measured with a higher sampling frequency, finite time effects could not be reduced, because time increments below the ME time scale must not be regarded in the KM analysis.
On time scales above the ME time scale, the process can almost perfectly be reconstructed by an OU process according to the classical overdamped Markov model of Brownian motion. The data quality does not allow to detect deviations from linearity in the drift term or a spatial dependence of the diffusion one could expect because the laser heats up the fluid.
The size of the measured diffusion coefficient is about 1.3 times smaller than the diffusion predicted by the Einstein-Stokes equation. For comparison, we also measured trajectories of different freely diffusing particles with our experimental setup. Averaged over all particles, the Einstein-Stokes equation was found to be valid. The fluctuations of diffusion constants among different particles, which can probably be traced back to fluctuations among particle radii, are large enough to explain the low diffusion found for the trapped particle.