[PYTHON] Algorithme de détection de pas de pas de Poisson

scénario

--Un algorithme qui détecte le point de pas de la trajectoire qui marche dans le processus de Poisson -Créé en référence à "Dynamique d'assemblage des microtubules à résolution moléculaire" Nature (2006) 709-712 ――Un style qui détecte d'abord le plus grand pas et ajoute progressivement de plus petits pas ―― La bonne chose à propos de cet algorithme est que les paramètres libres sont pratiquement nuls et qu'il n'y a pas besoin d'étalonnage.

• Testé à l'aide du poisson_stepper.py créé précédemment

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)

Résultat du calcul

• La ligne bleue est des données d'expérience simulées, la ligne noire convient ――Eh bien, il peut être détecté avec une bonne précision.