Create a poisson stepper with numpy.random
script
--Create a time series with Steps to test the Step detection algorithm --Add Gaussian noise to each point in the time series (conscious of single molecule experiment) --Pick up values according to probability distribution with numpy.random --Window average is easy if you use slices
poisson_stepper.py
#!/usr/bin/python
import numpy as np
import matplotlib.pyplot as plt
# Parameters
velocity = 60.0 # [nm/sec]
step_size = 50.0 # [nm]
k = velocity / step_size # decay constant for exponential distribution
dt = 1.0e-4 # [sec] time step for simulation
n_step = 10
sigma = 3.0 # [nm] amount of noise added to data
time_resolution = 0.01 # [sec]
n_data_averaged = int(time_resolution / dt)
# Produce trajectory
trajectory = []
for i in range(n_step+1):
trajectory += [i*step_size]*int(np.random.exponential() / dt)
# Average data
trajectory = [np.average(trajectory[i:i+n_data_averaged]) + np.random.normal(scale = sigma)
for i in range(0, len(trajectory), n_data_averaged)]
# Plot
plt.plot(trajectory)
plt.xlabel("time")
plt.ylabel("position")
plt.savefig("trajectory.png ")
Calculation result
Recommended Posts