[PYTHON] Poisson-Schrittschritt-Erkennungsalgorithmus

Skript

• Ein Algorithmus, der den Schrittpunkt der Trajektorie erkennt, der im Poisson-Prozess abläuft -Erstellt unter Bezugnahme auf "Assemblierungsdynamik von Mikrotubuli bei molekularer Auflösung" Nature (2006) 709-712 ――Ein Stil, der den größten Schritt zuerst erkennt und nach und nach kleinere Schritte hinzufügt ――Das Gute an diesem Algorithmus ist, dass die freien Parameter praktisch Null sind und keine Kalibrierung erforderlich ist.
• Getestet mit der zuvor erstellten poisson_stepper.py

chi2.py

#!/usr/bin/python                                                                                                                                                                                                                                                                          

import matplotlib.pyplot as plt
import numpy as np
import poisson_stepper
from scipy.optimize import curve_fit
import time
import timer


def detect_step(trajectory):
    # Initialize                                                                                                                                                                                                                                                                           
    min_chi2 = float("inf")
    min_t = 0
    delta_times_length = 0.0
    n = len(trajectory)

    for i in range(1, n - 1):
        # Array preparation                                                                                                                                                                                                                                                                
        left = trajectory[0:i]
        right = trajectory[i:n]

        # Calculate Chi2                                                                                                                                                                                                                                                                   
        current_chi2 = len(left) * np.var(left) + len(right) * np.var(right)

        # Update minimum values                                                                                                                                                                                                                                                            
        if min_chi2 > current_chi2:
            min_chi2 = current_chi2
            min_t = i
            delta_times_length = abs(np.average(left) - np.average(right)) * n

    return min_t, delta_times_length


def plot(trajectory, border, average_trajectory):

    plt.plot(trajectory)
    for i in range(1, len(border)):
        x = np.arange(border[i-1], border[i], 0.1)
        y_value = average_trajectory[i-1]
        y = [y_value for _ in range(0, len(x))]
        plt.plot(x, y, color='k', linewidth=2.0)

        if i == len(border) - 1:
            break

        y = np.arange(average_trajectory[i-1], average_trajectory[i], 0.1)
        x = [border[i] for _ in range(0, len(y))]
        plt.plot(x, y, color='k', linewidth=2.0)

    plt.show()


def detect_step_all_ranges(trajectory, border):
    max_delta_times_length = 0.0
    max_step_t = 0

    for iborder in range(0, len(border)-1):
	start = border[iborder] + 1
        end   = border[iborder+1] - 1
        trajectory_to_process = trajectory[start:end]
        step_t, delta_times_length = detect_step(trajectory_to_process)
        step_t += start

        if max_delta_times_length < delta_times_length:
            max_delta_times_length = delta_times_length
            max_step_t = step_t

    return max_step_t


def chi2(trajectory, t):

    # Calculate border                                                                                                                                                                                                                                                                     
    border = [0, len(trajectory)]

    for i in range(0, 10):
        border.append(detect_step_all_ranges(trajectory, border))
        border = sorted(border)

    # Calculate average trajectory                                                                                                                                                                                                                                                         
    average_trajectory = [np.average(trajectory[border[i]:border[i+1]])
                          for i in range(0, len(border) - 1)]

    # Plot                                                                                                                                                                                                                                                                                 
    plot(trajectory, border, average_trajectory)

    # Calculate step size                                                                                                                                                                                                                                                                  
    step_size = [average_trajectory[i+1] - average_trajectory[i]
                 for i in range(0, len(average_trajectory)-1)]

    return step_size


if __name__ == "__main__":
    trajectory, t = poisson_stepper.poisson_stepper(sigma=10.0)
    step_size = chi2(trajectory, t)

Berechnungsergebnis

• Blaue Linie sind simulierte Versuchsdaten, schwarze Linie passt ――Nun, es kann mit guter Genauigkeit erkannt werden.