Daily Archives: 2025年11月14日

Interpreting Binding Events in Single-Molecule Fluorescence Tracking: Conceptual Clarification and Practical Proposal

Leave a reply
  1. Concept: is the start of a track a binding event? From your plot: Each track appears suddenly at some time, then decays / fluctuates and disappears. The accumulated trace shows clear upward steps at those times. In typical kymograph / spot-tracking data: A track is created when the tracking algorithm first detects a fluorescent spot above background. That is usually when a fluorophore binds (or comes into focus / into the ROI). The end of a track is when the spot disappears, often due to bleaching or unbinding. So, conceptually: Yes, in many experiments the first time point of a track is a natural candidate for the binding time of that molecule. But there are a few caveats you might want to mention to Vero: Some molecules may already be bound before the camera started recording – in that case the first frame is not the true binding time, just the first time you see it. Detection thresholds may delay the detected “start” by a few frames. If two tracks overlap in time, the accumulated signal will have multiple steps, and the HMM sees “binding” whenever the sum intensity increases. That also explains your observations: Per-track HMM: each trace is relatively short and noisy → few clear transitions, so the HMM doesn’t detect many events. HMM on the accumulated signal: all tracks summed → big, clean steps when a new molecule appears → HMM picks up almost all track starts as binding events. So it’s very reasonable to treat track starts as a kind of “ground truth” binding time, at least for benchmarking.
  2. Improved text you can send to Vero Here’s a cleaned-up version of what you wrote, phrased as a question: From the plots of the individual tracks and the accumulated photon counts (see attached figure), I have a conceptual question about how to define the binding time. In the track-intensity file, each track appears at some time point and then decays. Would you consider the first time point of each track to be the “binding time” of that molecule? To me this looks like a very natural definition, because the track starts exactly when the signal first appears. Note: – When I run the HMM on individual tracks, only very few binding events are detected. – When I run the HMM on the accumulated signal, I can detect almost all of the track starts as binding events. This makes me think that the track starting times might be a good gold standard for binding time points. Does this interpretation match your understanding of the experiment? If you want I can also help you write a follow-up mail where you propose to use “track start = binding” and “track end = bleaching/unbinding” as ground truth for validating your HMM.

Here is a concise bilingual summary of your text explaining the idea and a cleaned-up proposal to send to Vero:

English Summary

  1. Concept: Is the start of a track a binding event?
  • Each fluorescent track appears suddenly, fluctuates, then disappears; cumulative traces show clear stepwise increments.
  • A track starts when a fluorophore binds or enters detection (fluorescent) focus, and ends due to bleaching or unbinding.
  • Thus, the start of a track commonly represents the molecule’s binding time.
  • Caveats:
    • Some molecules may be bound before recording starts, so first observed frame is not true binding time.
    • Detection thresholds may delay start detection.
    • Overlapping tracks cause summed signal steps, and HMM detects any increase as binding.
  • Per-track HMM sees few transitions due to noise and short traces; summed HMM on total signal detects most binding events.
  1. Cleaned-up question to Vero: “From individual and cumulative photon count plots (see attached), does the first time point of each track represent the molecule’s binding time?
  • HMM on single tracks detects very few events;
  • HMM on accumulated trace detects nearly all track starts as bindings;
  • This suggests track start time is a good benchmark for binding. Does this match your understanding? I can help draft a follow-up email proposing to use ‘track start = binding’ and ‘track end = bleaching/unbinding’ as ground truth for HMM validation.”

中文总结

  1. 概念解释：轨迹开始是否代表结合事件？
  • 每条轨迹突然出现，随后衰减或波动后消失，累积轨迹中出现明显的阶梯式上升。
  • 轨迹的开始通常是荧光标记分子结合或进入焦点检测区域的时刻，结束则因漂白或解离。
  • 因此轨迹开始时间通常被视为该分子结合的时间点。
  • 注意事项：
    • 某些分子可能在录像开始前已结合，首帧不是准确结合时间；
    • 检测阈值可能导致结合时间存在几十帧延迟；
    • 轨迹重叠造成信号累积，HMM会将信号上升视为结合事件。
  • 针对单条轨迹的HMM因噪声和轨迹短小，事件检测少；对累积信号的HMM检测到几乎所有结合。
  1. 整理后可以发给Vero的文本： “根据个别轨迹和累积光子计数图（见附件），我有个问题：是否可以将每条轨迹的首个时间点视为该分子的结合时刻？
  • 针对单轨迹的HMM只检测到少量结合事件；
  • 针对累积轨迹的HMM则几乎检测到所有轨迹期初的结合；
  • 这可能说明轨迹起始时间是一个可靠的结合时间‘黄金标准’。 这种理解符合你的看法吗？我也可以帮你起草一封邮件，提议用“轨迹开始=结合”，“轨迹结束=漂白/解离”作为验证HMM的真实标准。”

12345678

⁂

  1. https://pmc.ncbi.nlm.nih.gov/articles/PMC9150788/ ↩

  2. https://www.sciencedirect.com/science/article/pii/S0022283625003742 ↩

  3. https://phys.libretexts.org/Courses/University_of_California_Davis/Biophysics_241:_Membrane_Biology/05:_Experimental_Characterization_-_Spectroscopy_and_Microscopy/5.07:_Single_Molecule_Tracking ↩

  4. https://www.nature.com/articles/s43586-021-00038-x ↩

  5. https://elifesciences.org/articles/93183 ↩

  6. https://pmc.ncbi.nlm.nih.gov/articles/PMC5028128/ ↩

  7. https://academic.oup.com/nar/article/49/19/e112/6355882 ↩

  8. https://royalsocietypublishing.org/doi/10.1098/rsob.210383 ↩

使用 ICON HMM 在 Octave 中自动检测单分子 binding/bleach 事件并估计 dwell time

Leave a reply 
% detect_binding_bleach.m
% 用法（在终端）:
%   octave detect_binding_bleach.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv

%--------------------------------------------
% 基本设置
%--------------------------------------------
arg_list = argv();
if numel(arg_list) < 1
    error("Usage: octave detect_binding_bleach.m 
<track_intensity_csv>");
end
input_file = arg_list{1};

% 加载 HMM 采样器
%addpath("sampler_SRC");

% 如果安装了 statistics 包，用于 kmeans
try
    pkg load statistics;
catch
    warning("Could not load 'statistics' package. Make sure it is installed if kmeans is missing.");
end

%--------------------------------------------
% 读入 track intensity 文件（分号分隔）
%--------------------------------------------
fprintf("Reading file: %s\n", input_file);

fid = fopen(input_file, "r");
if fid < 0
    error("Cannot open file: %s", input_file);
end

% 第一行是注释头 "# track index;time (seconds);track intensity (photon counts)"
header_line = fgetl(fid);  % 忽略内容，只是读掉这一行

% 后面每行: track_index;time_sec;intensity
data = textscan(fid, "%f%f%f", "Delimiter", ";");
fclose(fid);

track_idx = data{1};
time_sec  = data{2};
counts    = data{3};

% 排序（确保同一个 track 内按时间排序）
[~, order] = sortrows([track_idx, time_sec], [1, 2]);
track_idx = track_idx(order);
time_sec  = time_sec(order);
counts    = counts(order);

tracks = unique(track_idx);
n_tracks = numel(tracks);

fprintf("Found %d tracks.\n", n_tracks);

%--------------------------------------------
% 结果结构体
%--------------------------------------------
results = struct( ...
    "track_id", {}, ...
    "binding_indices", {}, ...
    "binding_times", {}, ...
    "bleach_indices", {}, ...
    "bleach_times", {} );

%--------------------------------------------
% 循环每个 track，跑 HMM + 找事件
%--------------------------------------------
for ti = 1:n_tracks
    tr = tracks(ti);
    fprintf("\n===== Track %d =====\n", tr);

    sel = (track_idx == tr);
    t = time_sec(sel);
    z = counts(sel);
    z = z(:);      % 列向量

    %-------------------------
    % 1) 设定 ICON HMM 的参数
    %-------------------------
    opts = struct();

    % 超参数（和原 MATLAB 代码一致）
    opts.a = 1;
    opts.g = 1;

    opts.Q(1) = mean(z);
    opts.Q(2) = 1 / (std(z)^2 + eps);   % 防止除零
    opts.Q(3) = 0.1;
    opts.Q(4) = 0.00001;

    opts.M = 10;

    % 采样参数
    opts.dr_sk  = 1;
    opts.K_init = 50;

    flag_stat = true;
    flag_anim = false;    % 建议在 Octave 里关掉动画更稳定

    R = 1000;             % 采样次数，可根据需要调小/调大

    %-------------------------
    % 2) 运行采样器
    %-------------------------
    fprintf("  Running HMM sampler...\n");
    chain = chainer_main(z, R, [], opts, flag_stat, flag_anim);

    %-------------------------
    % 3) 做后验分析，得到平滑后的光子水平轨迹
    %-------------------------
    fr = 0.25;        % burn-in 比例
    dr = 2;           % sample 步距

    m_min = 0;
    m_max = 1;
    m_num = 25;

    [m_mod, m_red] = chainer_analyze_means(chain, fr, dr, m_min, m_max, m_num, z);

    % m_mod 是一个和 z 同长度的向量（通常是 [0,1] 范围的归一化水平）
    m_mod = m_mod(:);

    %-------------------------
    % 4) 把 m_mod 聚类成 3 个光子强度状态
    %    状态 1 = 最低强度 (dark / unbound / bleached)
    %    状态 3 = 最高强度 (bound)
    %-------------------------
    K = 3;   % 你可以改成 2 或其他
    try
        [idx_raw, centers] = kmeans(m_mod, K);
    catch
        % 如果没 kmeans，就简单按数值分位数来分三类
        warning("kmeans not available, using simple quantile-based clustering.");
        q1 = quantile(m_mod, 1/3);
        q2 = quantile(m_mod, 2/3);
        idx_raw = ones(size(m_mod));
        idx_raw(m_mod > q1 & m_mod <= q2) = 2;
        idx_raw(m_mod > q2) = 3;
        centers = zeros(K,1);
        for kk = 1:K
            centers(kk) = mean(m_mod(idx_raw == kk));
        end
    end

    % 根据中心值从小到大重排状态编号
    [~, order_centers] = sort(centers);
    state_seq = zeros(size(idx_raw));
    for k = 1:K
        state_seq(idx_raw == order_centers(k)) = k;
    end

    low_state  = 1;
    high_state = K;

    %-------------------------
    % 5) 检测 binding / bleach 跳变
    %    binding:  low -> high
    %    bleach : high -> low
    %-------------------------
    s = state_seq(:);
    % 前一时刻和后一时刻
    s_prev = s(1:end-1);
    s_next = s(2:end);

    bind_idx   = find(s_prev == low_state  & s_next == high_state) + 1;
    bleach_idx = find(s_prev == high_state & s_next == low_state) + 1;

    bind_times   = t(bind_idx);
    bleach_times = t(bleach_idx);

    fprintf("  Found %d binding event(s) and %d bleaching event(s).\n", ...
            numel(bind_idx), numel(bleach_idx));

    % 存入结果
    results(ti).track_id        = tr;
    results(ti).binding_indices = bind_idx;
    results(ti).binding_times   = bind_times;
    results(ti).bleach_indices  = bleach_idx;
    results(ti).bleach_times    = bleach_times;
end

%--------------------------------------------
% 6) 把结果写成 CSV
%--------------------------------------------
out_csv = "binding_bleach_events.csv";
fid_out = fopen(out_csv, "w");
if fid_out < 0
    error("Cannot open output file: %s", out_csv);
end

fprintf(fid_out, "track_index,event_type,sample_index,time_seconds\n");
for ti = 1:numel(results)
    tr = results(ti).track_id;

    % binding
    for k = 1:numel(results(ti).binding_indices)
        fprintf(fid_out, "%d,binding,%d,%.6f\n", tr, ...
                results(ti).binding_indices(k), ...
                results(ti).binding_times(k));
    end

    % bleach
    for k = 1:numel(results(ti).bleach_indices)
        fprintf(fid_out, "%d,bleach,%d,%.6f\n", tr, ...
                results(ti).bleach_indices(k), ...
                results(ti).bleach_times(k));
    end
end
fclose(fid_out);

fprintf("\nAll done. Events written to: %s\n", out_csv);
  • 运行：
octave detect_binding_bleach.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv
  1. 保留原始 ICON 后验分析：
    • chainer_export 导出样本
    • chainer_analyze_means 得到 m_mod
    • chainer_analyze_transitions 得到 m_edges, p_mean, d_dist
    • chainer_analyze_drift 得到 y_mean, y_std
  2. 使用 binding/bleach 时间点 计算 dwell time 分布并做单指数 rate 拟合，并直接画直方图 + 指数拟合曲线。
  3. 输出一个 CSV 表，包含所有 binding/bleach 时间点（这是最重要的输出）。

说明：

  • 代码仍然基于你已有的 detect_binding_bleach.m 结构，假定当前目录下有 sampler_SRC（含 chainer_main.m、chainer_analyze_means.m、chainer_analyze_transitions.m、chainer_analyze_drift.m 等）。910
  • Dwell time 这里使用“连续 binding 事件之间的时间差”来估计 bound-state dwell time（简化版）；如果 bleaching 明确是轨迹末端终止，也可以用 bleach_time - last_binding_time 作为一个 right‑censored dwell（这里暂不做 censoring 修正，只输出原始分布和单指数拟合）。

完整脚本：detect_binding_bleach_dwell_simple.m

% detect_binding_bleach_dwell.m
%
% 用法（终端）:
%   octave detect_binding_bleach_dwell.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv
%
% 输入 CSV 格式（分号分隔）:
%   # track index;time (seconds);track intensity (photon counts)
%   1;0.00;123
%   1;0.05;118
%   2;0.00; 95
%   ...

%--------------------------------------------
% 0. 基本设置
%--------------------------------------------
arg_list = argv();

if numel(arg_list) < 1
    error("Usage: octave detect_binding_bleach_dwell.m 
<input_csv>");
end

input_file = arg_list{1};

% 加载 HMM 采样器代码（确保 sampler_SRC 在当前目录下）
addpath("sampler_SRC");

% 如果安装了 statistics 包，用于 kmeans
try
    pkg load statistics;
catch
    warning("Could not load 'statistics' package. Make sure it is installed if kmeans is missing.");
end

%--------------------------------------------
% 1. 读入 track intensity 文件（分号分隔）
%--------------------------------------------
fprintf("Reading file: %s\n", input_file);

fid = fopen(input_file, "r");
if fid < 0
    error("Cannot open file: %s", input_file);
end

% 第一行是注释头 "# track index;time (seconds);track intensity (photon counts)"
header_line = fgetl(fid); % 忽略内容，只是读掉这一行

% 后面每行: track_index;time_sec;intensity
data = textscan(fid, "%f%f%f", "Delimiter", ";");
fclose(fid);

track_idx = data{1};
time_sec  = data{2};
counts    = data{3};

% 排序（确保同一个 track 内按时间排序）
[~, order] = sortrows([track_idx, time_sec], [1, 2]);
track_idx = track_idx(order);
time_sec  = time_sec(order);
counts    = counts(order);

tracks   = unique(track_idx);
n_tracks = numel(tracks);
fprintf("Found %d tracks.\n", n_tracks);

%--------------------------------------------
% 2. 结果结构体（binding / bleach 时间点）
%--------------------------------------------
results = struct( ...
    "track_id",        {}, ...
    "binding_indices", {}, ...
    "binding_times",   {}, ...
    "bleach_indices",  {}, ...
    "bleach_times",    {} );

% 用于 dwell time 收集的数组
all_binding_times = [];   % 所有 track 的 binding time 合并
all_bleach_times  = [];   % 所有 track 的 bleach time 合并

%--------------------------------------------
% 3. 循环每个 track，跑 HMM + ICON 分析 + 事件检测
%--------------------------------------------
for ti = 1:n_tracks

    tr = tracks(ti);
    fprintf("\n===== Track %d =====\n", tr);

    sel = (track_idx == tr);
    t   = time_sec(sel);
    z   = counts(sel);
    z   = z(:); % 列向量

    %-------------------------
    % 3.1 设定 ICON HMM 的参数
    %-------------------------
    opts        = struct();
    % 超参数（与 level_analysis.m 一致）
    opts.a      = 1;
    opts.g      = 1;
    opts.Q(1)   = mean(z);
    opts.Q(2)   = 1 / (std(z)^2 + eps); % 防止 std(z)=0 除零
    opts.Q(3)   = 0.1;
    opts.Q(4)   = 0.00001;
    opts.M      = 10;

    % 采样参数
    opts.dr_sk  = 1;
    opts.K_init = 50;

    flag_stat   = true;
    flag_anim   = false;  % 建议在 Octave 里关掉动画
    R           = 1000;   % 采样次数，可按需求调整

    %-------------------------
    % 3.2 运行采样器（ICON HMM）
    %-------------------------
    fprintf(" Running HMM sampler...\n");
    chain = chainer_main(z, R, [], opts, flag_stat, flag_anim);

    %-------------------------
    % 3.3 ICON 后验分析（与 level_analysis.m 等价）
    %-------------------------
    fr    = 0.25; % burn-in 比例
    dr    = 2;    % sample 步距
    m_min = 0;
    m_max = 1;
    m_num = 25;

    % (1) 均值轨迹：m_mod
    [m_mod, m_red] = chainer_analyze_means(chain, fr, dr, m_min, m_max, m_num, z);
    m_mod = m_mod(:);

    % (2) 导出 samples（为每条 track 单独存一个文件）
    sample_file = sprintf("samples_track_%d", tr);
    chainer_export(chain, fr, dr, sample_file, "mat");

    % (3) 转移概率 / 跃迁统计
    [m_edges, p_mean, d_dist] = chainer_analyze_transitions( ...
        chain, fr, dr, m_min, m_max, m_num, true);

    % (4) 漂移分析
    [y_mean, y_std] = chainer_analyze_drift(chain, fr, dr, z);

    % 你可以根据需要保存这些 ICON 分析结果
    mat_out = sprintf("icon_analysis_track_%d.mat", tr);
    save(mat_out, "m_mod", "m_red", "m_edges", "p_mean", "d_dist", ...
                  "y_mean", "y_std", "t", "z");

    %-------------------------
    % 3.4 把 m_mod 聚类成 3 个光子强度状态（用于事件检测）
    %-------------------------
    K = 3; % 可按物理需求修改

    try
        [idx_raw, centers] = kmeans(m_mod, K);
    catch
        warning("kmeans not available, using simple quantile-based clustering.");
        q1 = quantile(m_mod, 1/3);
        q2 = quantile(m_mod, 2/3);

        idx_raw = ones(size(m_mod));
        idx_raw(m_mod > q1 & m_mod <= q2) = 2;
        idx_raw(m_mod > q2)              = 3;

        centers = zeros(K,1);
        for kk = 1:K
            centers(kk) = mean(m_mod(idx_raw == kk));
        end
    end

    % 根据中心值从小到大重排状态编号，使 state_seq = 1..K 对应 low->high
    [~, order_centers] = sort(centers);
    state_seq = zeros(size(idx_raw));
    for k = 1:K
        state_seq(idx_raw == order_centers(k)) = k;
    end

    low_state  = 1;
    high_state = K;

    %-------------------------
    % 3.5 检测 binding / bleaching
    %
    %   binding: low_state -> high_state
    %   bleach : high_state -> low_state
    %-------------------------
    s      = state_seq(:);
    s_prev = s(1:end-1);
    s_next = s(2:end);

    bind_idx   = find(s_prev == low_state  & s_next == high_state) + 1;
    bleach_idx = find(s_prev == high_state & s_next == low_state) + 1;

    bind_times   = t(bind_idx);
    bleach_times = t(bleach_idx);

    fprintf(" Found %d binding event(s) and %d bleaching event(s).\n", ...
            numel(bind_idx), numel(bleach_idx));

    % 存入结果结构体
    results(ti).track_id        = tr;
    results(ti).binding_indices = bind_idx;
    results(ti).binding_times   = bind_times;
    results(ti).bleach_indices  = bleach_idx;
    results(ti).bleach_times    = bleach_times;

    % 用于全局 dwell-time 统计
    all_binding_times = [all_binding_times; bind_times(:)];
    all_bleach_times  = [all_bleach_times;  bleach_times(:)];

end

%--------------------------------------------
% 4. 输出 binding/bleach 时间点的总表（最重要输出）
%--------------------------------------------
out_csv = "binding_bleach_events.csv";
fid_out = fopen(out_csv, "w");
if fid_out < 0
    error("Cannot open output file: %s", out_csv);
end

fprintf(fid_out, "track_index,event_type,sample_index,time_seconds\n");

for ti = 1:numel(results)
    tr = results(ti).track_id;

    % binding
    for k = 1:numel(results(ti).binding_indices)
        fprintf(fid_out, "%d,binding,%d,%.6f\n", tr, ...
                results(ti).binding_indices(k), ...
                results(ti).binding_times(k));
    end

    % bleach
    for k = 1:numel(results(ti).bleach_indices)
        fprintf(fid_out, "%d,bleach,%d,%.6f\n", tr, ...
                results(ti).bleach_indices(k), ...
                results(ti).bleach_times(k));
    end
end

fclose(fid_out);
fprintf("\n[INFO] Binding/bleach events written to: %s\n", out_csv);

%--------------------------------------------
% 5. Dwell-time 直方图 + rate 拟合
%
%   简化假设：
%   - 把每个 track 中连续 binding 事件之间的时间差视为
%     bound state dwell time（中间没有 bleaching，就重新 binding）。
%   - 如果有明确的最后一个 binding 到 bleaching 的时间，
%     你也可以扩展这里把它们当作 dwell time 加进去。
%--------------------------------------------

% 5.1 收集 bound-state dwell times（简单用 binding 时间差）
dwell_times = [];

for ti = 1:numel(results)
    bt = sort(results(ti).binding_times(:)); % 单个 track 的 binding times
    if numel(bt) >= 2
        dt = diff(bt);    % 相邻 binding 之间的时间差
        dwell_times = [dwell_times; dt];
    end

    % 可选：如果你想把最后一个 binding -> 第一个 bleach 也算入 dwell:
    % if ~isempty(results(ti).binding_times) && ~isempty(results(ti).bleach_times)
    %     last_binding  = max(results(ti).binding_times);
    %     first_bleach  = min(results(ti).bleach_times);
    %     if first_bleach > last_binding
    %         dwell_times = [dwell_times; first_bleach - last_binding];
    %     end
    % end
end

if isempty(dwell_times)
    fprintf("[WARN] No sufficient binding events to compute dwell-time distribution.\n");
else
    % 5.2 绘制 dwell-time 直方图
    figure;
    hold on;
    nbins = max(10, round(sqrt(numel(dwell_times)))); % 简单经验 bin 数
    [counts_hist, edges] = histcounts(dwell_times, nbins);
    centers = (edges(1:end-1) + edges(2:end)) / 2;

    bar(centers, counts_hist, "hist");
    xlabel("Dwell time (s)");
    ylabel("Counts");
    title("Bound-state dwell-time histogram");

    % 5.3 拟合单指数: p(t) ~ (1/tau) * exp(-t/tau)
    % 使用简单的最小二乘拟合 log(count) 对 t
    valid = counts_hist > 0;
    t_fit = centers(valid);
    y_fit = log(counts_hist(valid));

    % 线性拟合: log(N) = a + b * t, 其中 b ~= -1/tau
    p = polyfit(t_fit, y_fit, 1);
    a = p(1);
    b = p(2);

    tau_est = -1 / b;
    fprintf("[INFO] Fitted single-exponential dwell time tau = %.4f s\n", tau_est);

    % 画出拟合曲线（归一到 histogram 的总 counts 约束）
    t_plot = linspace(min(dwell_times), max(dwell_times), 200);
    % 预测密度 ~ exp(a + b t)，但我们要匹配 histogram 的离散计数
    y_model = exp(a + b * t_plot);

    plot(t_plot, y_model, "r-", "LineWidth", 2);
    legend("Histogram", "Single-exp fit");
    hold off;

    % 可选：保存图像
    print("-dpng", "dwell_time_histogram_fit.png");
    fprintf("[INFO] Dwell-time histogram and fit saved to dwell_time_histogram_fit.png\n");
end

fprintf("\nAll done.\n");

使用与扩展建议

  • 运行：
octave detect_binding_bleach_dwell_simple.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv

会生成：

  • binding_bleach_events.csv：所有 track 的 binding/bleach 时间点（你最关心的表格）。
  • icon_analysis_track_X.mat：每个 track 的 ICON 分析结果（m_mod, m_edges, p_mean, d_dist, y_mean, y_std 等）。
  • 如果有足够事件，生成 dwell_time_histogram_fit.png，包含 dwell-time 直方图和单指数拟合曲线。
    • dwell time 定义： 当前用的是“相邻 binding 时间点之差”作为 dwell；你如果有更明确的物理模型（例如每条轨迹只有 1 次 binding + 1 次 bleach），可以把注释里那段 last_binding -> first_bleach 的逻辑打开或进一步改写成更严格的两态 dwell 定义。
    • rate constant： 对于单指数 $p(t) \propto \exp(-t/\tau)$，tau_est 是平均 dwell time，解离 rate 大致是 $k{\text{off}} = 1/\tau{\text{est}}$。如果你需要可以在脚本里再打印一条 k_off_est。

如果你希望 dwell-time 统计严格区分 “bound state dwell vs unbound state dwell”，或者想直接从 p_mean / d_dist 中做 rate 矩阵估计，可以在这份脚本基础上再往前一步；可以一起设计一个更“教科书版”的 HMM‑rate 推断流程。

特性概览：

  • 使用 ICON HMM（chainer_main）分析每个 track。9
  • 做原始 level_analysis.m 的四类后验分析：chainer_export、chainer_analyze_means、chainer_analyze_transitions、chainer_analyze_drift。9
  • 基于 m_mod 聚类成 3 个强度状态，检测 binding / bleaching 事件。10
  • 输出 binding_bleach_events.csv，包含所有 binding/bleach 时间点。10
  • 用 binding 时间点计算 dwell-time 分布，做单指数拟合；绘图部分只在内存中生成，用 print + close 保存 PNG（类似 R 的 dev.off()），并使用 hist 而不是 histcounts。10

将以下内容保存为 detect_binding_bleach_dwell.m 即可直接运行。

detect_binding_bleach_dwell.m 完整代码

% detect_binding_bleach_dwell.m
%
% 用法（终端）:
%   octave detect_binding_bleach_dwell.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv
%
% 输入 CSV 格式（分号分隔）:
%   # track index;time (seconds);track intensity (photon counts)
%   1;0.00;123
%   1;0.05;118
%   2;0.00; 95
%   ...

%--------------------------------------------
% 0. 基本设置
%--------------------------------------------
arg_list = argv();

if numel(arg_list) < 1
    error("Usage: octave detect_binding_bleach_dwell.m 
<input_csv>");
end

input_file = arg_list{1};

% 加载 HMM 采样器代码（确保 sampler_SRC 在当前目录下）
addpath("sampler_SRC");

% 如果安装了 statistics 包，用于 kmeans
try
    pkg load statistics;
catch
    warning("Could not load 'statistics' package. Make sure it is installed if kmeans is missing.");
end

%--------------------------------------------
% 1. 读入 track intensity 文件（分号分隔）
%--------------------------------------------
fprintf("Reading file: %s\n", input_file);

fid = fopen(input_file, "r");
if fid < 0
    error("Cannot open file: %s", input_file);
end

% 第一行是注释头 "# track index;time (seconds);track intensity (photon counts)"
header_line = fgetl(fid); % 忽略内容，只是读掉这一行

% 后面每行: track_index;time_sec;intensity
data = textscan(fid, "%f%f%f", "Delimiter", ";");
fclose(fid);

track_idx = data{1};
time_sec  = data{2};
counts    = data{3};

% 排序（确保同一个 track 内按时间排序）
[~, order] = sortrows([track_idx, time_sec], [1, 2]);
track_idx = track_idx(order);
time_sec  = time_sec(order);
counts    = counts(order);

tracks   = unique(track_idx);
n_tracks = numel(tracks);
fprintf("Found %d tracks.\n", n_tracks);

%--------------------------------------------
% 2. 结果结构体（binding / bleach 时间点）
%--------------------------------------------
results = struct( ...
    "track_id",        {}, ...
    "binding_indices", {}, ...
    "binding_times",   {}, ...
    "bleach_indices",  {}, ...
    "bleach_times",    {} );

% 用于 dwell time 收集的数组（跨所有 track）
all_binding_times = [];   % 所有 track 的 binding time 合并
all_bleach_times  = [];   % 所有 track 的 bleach time 合并

%--------------------------------------------
% 3. 循环每个 track，跑 HMM + ICON 分析 + 事件检测
%--------------------------------------------
for ti = 1:n_tracks

    tr = tracks(ti);
    fprintf("\n===== Track %d =====\n", tr);

    sel = (track_idx == tr);
    t   = time_sec(sel);
    z   = counts(sel);
    z   = z(:); % 列向量

    %-------------------------
    % 3.1 设定 ICON HMM 的参数
    %-------------------------
    opts        = struct();
    % 超参数（与 level_analysis.m 一致）
    opts.a      = 1;
    opts.g      = 1;
    opts.Q(1)   = mean(z);
    opts.Q(2)   = 1 / (std(z)^2 + eps); % 防止 std(z)=0 除零
    opts.Q(3)   = 0.1;
    opts.Q(4)   = 0.00001;
    opts.M      = 10;

    % 采样参数
    opts.dr_sk  = 1;
    opts.K_init = 50;

    flag_stat   = true;
    flag_anim   = false;  % 建议在 Octave 里关掉动画
    R           = 1000;   % 采样次数，可按需求调整

    %-------------------------
    % 3.2 运行采样器（ICON HMM）
    %-------------------------
    fprintf(" Running HMM sampler...\n");
    chain = chainer_main(z, R, [], opts, flag_stat, flag_anim);

    %-------------------------
    % 3.3 ICON 后验分析（与 level_analysis.m 等价）
    %-------------------------
    fr    = 0.25; % burn-in 比例
    dr    = 2;    % sample 步距
    m_min = 0;
    m_max = 1;
    m_num = 25;

    % (1) 均值轨迹：m_mod
    [m_mod, m_red] = chainer_analyze_means(chain, fr, dr, m_min, m_max, m_num, z);
    m_mod = m_mod(:);

    % (2) 导出 samples（为每条 track 单独存一个文件）
    sample_file = sprintf("samples_track_%d", tr);
    chainer_export(chain, fr, dr, sample_file, "mat");
    fprintf("%s --- Exported\n", [sample_file ".mat"]);

    % (3) 转移概率 / 跃迁统计
    [m_edges, p_mean, d_dist] = chainer_analyze_transitions( ...
        chain, fr, dr, m_min, m_max, m_num, true);

    % (4) 漂移分析
    [y_mean, y_std] = chainer_analyze_drift(chain, fr, dr, z);

    % 保存这些 ICON 分析结果（可选）
    mat_out = sprintf("icon_analysis_track_%d.mat", tr);
    save(mat_out, "m_mod", "m_red", "m_edges", "p_mean", "d_dist", ...
                  "y_mean", "y_std", "t", "z");

    %-------------------------
    % 3.4 把 m_mod 聚类成 3 个光子强度状态（用于事件检测）
    %-------------------------
    K = 3; % 可按物理需求修改

    try
        [idx_raw, centers] = kmeans(m_mod, K);
    catch
        warning("kmeans not available, using simple quantile-based clustering.");
        q1 = quantile(m_mod, 1/3);
        q2 = quantile(m_mod, 2/3);

        idx_raw = ones(size(m_mod));
        idx_raw(m_mod > q1 & m_mod <= q2) = 2;
        idx_raw(m_mod > q2)              = 3;

        centers = zeros(K,1);
        for kk = 1:K
            centers(kk) = mean(m_mod(idx_raw == kk));
        end
    end

    % 根据中心值从小到大重排状态编号，使 state_seq = 1..K 对应 low->high
    [~, order_centers] = sort(centers);
    state_seq = zeros(size(idx_raw));
    for k = 1:K
        state_seq(idx_raw == order_centers(k)) = k;
    end

    low_state  = 1;
    high_state = K;

    %-------------------------
    % 3.5 检测 binding / bleaching
    %
    %   binding: low_state -> high_state
    %   bleach : high_state -> low_state
    %-------------------------
    s      = state_seq(:);
    s_prev = s(1:end-1);
    s_next = s(2:end);

    bind_idx   = find(s_prev == low_state  & s_next == high_state) + 1;
    bleach_idx = find(s_prev == high_state & s_next == low_state) + 1;

    bind_times   = t(bind_idx);
    bleach_times = t(bleach_idx);

    fprintf("Found %d binding event(s) and %d bleaching event(s).\n", ...
            numel(bind_idx), numel(bleach_idx));

    % 存入结果结构体
    results(ti).track_id        = tr;
    results(ti).binding_indices = bind_idx;
    results(ti).binding_times   = bind_times;
    results(ti).bleach_indices  = bleach_idx;
    results(ti).bleach_times    = bleach_times;

    % 用于全局 dwell-time 统计
    all_binding_times = [all_binding_times; bind_times(:)];
    all_bleach_times  = [all_bleach_times;  bleach_times(:)];

end

%--------------------------------------------
% 4. 输出 binding/bleach 时间点的总表（最重要输出）
%--------------------------------------------
out_csv = "binding_bleach_events.csv";
fid_out = fopen(out_csv, "w");
if fid_out < 0
    error("Cannot open output file: %s", out_csv);
end

fprintf(fid_out, "track_index,event_type,sample_index,time_seconds\n");

for ti = 1:numel(results)
    tr = results(ti).track_id;

    % binding
    for k = 1:numel(results(ti).binding_indices)
        fprintf(fid_out, "%d,binding,%d,%.6f\n", tr, ...
                results(ti).binding_indices(k), ...
                results(ti).binding_times(k));
    end

    % bleach
    for k = 1:numel(results(ti).bleach_indices)
        fprintf(fid_out, "%d,bleach,%d,%.6f\n", tr, ...
                results(ti).bleach_indices(k), ...
                results(ti).bleach_times(k));
    end
end

fclose(fid_out);
fprintf("\n[INFO] Binding/bleach events written to: %s\n", out_csv);

%--------------------------------------------
% 5. Dwell-time 直方图 + rate 拟合
%
%   简化假设：
%   - 把每个 track 中连续 binding 事件之间的时间差视为
%     bound state dwell time（中间没有 bleaching，就重新 binding）。
%   - 如果有明确的最后一个 binding 到 bleaching 的时间，
%     可以扩展此处加入（当前先保留简单版本）。
%--------------------------------------------

% 5.1 收集 bound-state dwell times（简单用 binding 时间差）
dwell_times = [];

for ti = 1:numel(results)
    bt = sort(results(ti).binding_times(:)); % 单个 track 的 binding times
    if numel(bt) >= 2
        dt = diff(bt);    % 相邻 binding 之间的时间差
        dwell_times = [dwell_times; dt];
    end

    % 可选：如果你想把最后一个 binding -> 第一个 bleach 也算入 dwell:
    % if ~isempty(results(ti).binding_times) && ~isempty(results(ti).bleach_times)
    %     last_binding  = max(results(ti).binding_times);
    %     first_bleach  = min(results(ti).bleach_times);
    %     if first_bleach > last_binding
    %         dwell_times = [dwell_times; first_bleach - last_binding];
    %     end
    % end
end

if isempty(dwell_times)
    fprintf("[WARN] No sufficient binding events to compute dwell-time distribution.\n");
else
    % 5.2 绘制 dwell-time 直方图（只在内存中画，不弹窗）
    fig = figure("visible", "off");
    hold on;

    nbins = max(10, round(sqrt(numel(dwell_times))));

    % 使用 hist（Octave 自带）代替 histcounts
    [counts_hist, centers] = hist(dwell_times, nbins);

    bar(centers, counts_hist, "hist");
    xlabel("Dwell time (s)");
    ylabel("Counts");
    title("Bound-state dwell-time histogram");

    % 5.3 拟合单指数: p(t) ~ (1/tau) * exp(-t/tau)
    valid = counts_hist > 0;
    t_fit = centers(valid);
    y_fit = log(counts_hist(valid));

    p = polyfit(t_fit, y_fit, 1);
    a = p(1);
    b = p(2);

    tau_est = -1 / b;
    fprintf("[INFO] Fitted single-exponential dwell time tau = %.4f s\n", tau_est);

    t_plot = linspace(min(dwell_times), max(dwell_times), 200);
    y_model = exp(a + b * t_plot);

    plot(t_plot, y_model, "r-", "LineWidth", 2);
    legend("Histogram", "Single-exp fit");
    hold off;

    % 保存图像并关闭 figure（类似 R 的 dev.off()）
    out_png = "dwell_time_histogram_fit.png";
    print(fig, out_png, "-dpng");
    close(fig);

    fprintf("[INFO] Dwell-time histogram and fit saved to %s\n", out_png);
end

fprintf("\nAll done.\n");

运行方式：

octave detect_binding_bleach_dwell.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv

如有下一步你想基于 d_dist 做更“纯 HMM”的 dwell 分布 / rate 矩阵，也可以在这个脚本上再加一节专门分析。

⁂

level_analysis.m vs detect_binding_bleach.m

这两个脚本在核心 HMM 采样算法上是一致的，但整体流程和用途并不完全相同。level_analysis.m只是做“ICON HMM 拟合 + 若干通用后验分析”，而 detect_binding_bleach.m在此基础上，额外做了“多 track 输入 + 聚类成离散状态 + 自动识别 binding/bleach 跳变 + 输出 CSV”。因此不能说是完全相同算法的简单移植，而是“在相同 HMM 内核上扩展出的一个更完整事件检测流程”。910

下面分两部分说明：先列出关键差异，然后用中文系统讲解 detect_binding_bleach.m 的整个处理流程（包括和原始 level_analysis.m 共享的 HMM 部分）。109

一、算法差异概览

  1. 输入数据形式不同
    • level_analysis.m：直接从 photoncount.txt 读入一个向量 z，视为单一时间序列（一个 DNA tether 的 photon counts）。9
    • detect_binding_bleach.m：从命令行传入一个分号分隔的 CSV，结构为 track_index;time;intensity，先按 (track_index, time) 排序，然后对每个 track 单独取出 z 与时间 t 做分析，相当于一次运行中批量处理多条轨迹。10
  2. HMM 部分是否一致
    • 两者都使用 chainer_main 这一 ICON HMM 采样器，参数设置非常接近：
      • opts.a = 1; opts.g = 1;（浓度参数）。910
      • opts.Q(1) = mean(z);（均值超分布的均值）。109
      • opts.Q(2) = 1/std(z)^2; vs 1 / (std(z)^2 + eps)（只是在 Octave 版本里加了 eps 防除零，对数值稳定性有轻微改进，但理论上等价）。910
      • opts.Q(3) = 0.1; opts.Q(4) = 0.00001; opts.M = 10; 完全一致。109
      • 采样设置：opts.dr_sk = 1; opts.K_init = 50; R = 1000; 也一致。910
    • 因此，给定同一 z，ICON HMM 的抽样和后验均值轨迹的生成逻辑是相同的，只是 Octave 版本做了防数值异常的小修正。整体可以认为 HMM 核心算法是等价的。109
  3. 后验分析部分的差异
    • level_analysis.m：侧重“通用分析”，调用：
      • chainer_export 导出样本。9
      • chainer_analyze_means 得到状态均值轨迹 m_mod。9
      • chainer_analyze_transitions 得到转移概率矩阵、边界等。9
      • chainer_analyze_drift 得到漂移轨迹 y_mean, y_std。9
    • detect_binding_bleach.m：只调用 chainer_analyze_means 得到 m_mod，然后完全走另一条分析线：
      • 用 kmeans 或分位数把 m_mod 聚成 3 个强度状态（低、中、高）。10
      • 通过状态序列中的 low→high 和 high→low 跳变来定义 binding 和 bleaching 事件，并输出事件时间。10
    • 所以 HMM 后处理的目标和输出是完全不同的：一个是做统计分析、可视化；另一个是做自动事件检测。
  4. 输出结果不同
    • level_analysis.m：
      • 导出 MCMC 样本 samples.mat，方便后续任意分析；同时有若干变量在 workspace 中（m_mod, p_mean, d_dist, y_mean 等）。9
    • detect_binding_bleach.m：
      • 对每条 track 计算并存储 binding_indices/binding_times 和 bleach_indices/bleach_times，最后全部写入统一的 binding_bleach_events.csv，用于后续统计或可视化。10
  5. 附加逻辑（聚类与事件判定）是新加的
    • kmeans（或分位数替代）、状态重标号、状态序列差分、事件 CSV 输出等逻辑，在 level_analysis.m 中完全不存在，是在 Octave 版本里新增的功能层。也就是说，detect_binding_bleach.m 不只是移植，而是在 HMM 内核基础上的一个“面向 binding/bleach 识别”的上层算法。109

结论：

  • HMM 模型和采样部分与原 MATLAB 代码在理论上是相同的，只做了数值稳定性改动。
  • 整体脚本逻辑不相同，Octave 版多了一整套“轨迹遍历 + 状态聚类 + binding/bleach 检测 + CSV 输出”的算法层。

二、用中文说明 detect_binding_bleach.m 的完整流程（含与原脚本共用的 HMM 思路）

下面把 detect_binding_bleach.m 当作主角，用中文分步骤解释其处理流程，并在相关步骤指出它与原 level_analysis.m 的对应关系。这样方便你检查是否满足你对“算法一致性”的期待。109

1. 命令行输入与文件读取

  • 脚本通过 argv() 读取命令行参数，第一个参数是输入文件名，例如： octave detect_binding_bleach.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv。10
  • 输入文件是以分号分隔的文本，第一行为注释头，后面每行包含三列：
    • track_index：轨迹编号（一个 DNA tether 或一个粒子）。
    • time_seconds：该轨迹下某一帧的时间点。
    • track_intensity：该时间点的 photon counts。
  • 读入后，把三列分别存成 track_idx, time_sec, counts。然后用 sortrows([track_idx, time_sec]) 按 “轨迹号优先，再按时间” 排好顺序，保证每条轨迹都是时间有序的。之后提取所有不同的 tracks，准备对每一条轨迹单独做分析。10

这一步在原 level_analysis.m 中要简单得多，那里只是 z = load('photoncount.txt');，默认只处理单一系列，没有 track 结构。9

2. 对每条轨迹循环，构建输入序列 z

  • 外层 for ti = 1:n_tracks 遍历每一个 track。
  • 对于当前轨迹 tr，用 sel = (track_idx == tr) 挑出这一条轨迹对应的时间和强度，得到：
    • t = time_sec(sel)：该轨迹的时间轴。
    • z = counts(sel)：该轨迹的 photon count 序列。
  • 把 z 转成列向量，作为后续 HMM 的观测数据。

到这里为止，z 的角色与原脚本中的 z = load('photoncount.txt') 完全一致，只是这里是“多轨迹版本”的 z。910

3. 设置 ICON HMM 的超参数和采样参数

这一段是两份脚本的核心重合部分，也是你最关心的“算法一致性”所在。910

  • 先构建 opts 结构体，用于向 chainer_main 传参。
  • 超参数设置为：
    • opts.a = 1;：控制 HMM 中状态转换的 Dirichlet 浓度（类似“转移稀疏度”的先验强度）。
    • opts.g = 1;：控制状态先验分布的基本参数（base measure 的浓度）。
    • opts.Q(1) = mean(z);：假设状态均值的先验均值等于观测平均值。
    • opts.Q(2) = 1 / (std(z)^2 + eps);：状态均值先验的精度，与观测标准差成反比；加上 eps 是为了避免在极端情形下 std(z) = 0 导致除零。原 MATLAB 版本是 1/std(z)^2，数学上是一致的设定。
    • opts.Q(3) = 0.1;：方差逆数的 Gamma 先验的 shape。
    • opts.Q(4) = 0.00001;：方差逆数的 Gamma 先验的 scale。
    • opts.M = 10;：插值节点数（ICON 模型内部数值近似需要）。
  • 采样相关的参数：
    • opts.dr_sk = 1;：保存样本的间隔步长。
    • opts.K_init = 50;：初始假设的隐状态数（实际有效状态数由采样自动调整）。
    • R = 1000;：采样迭代次数。

这些设置在 level_analysis.m 中完全一致，只是缺少 eps。所以从贝叶斯 HMM 模型结构和先验设定来看，两者是同一个 ICON 模型。109

4. 运行 HMM 采样器，得到 MCMC 链

  • 调用 chain = chainer_main(z, R, [], opts, flag_stat, flag_anim); 来运行 ICON HMM 的 Gibbs sampler 或类似的 MCMC 算法。
    • z 是观测序列。
    • R 是采样轮数。
    • 第三个参数为空 []，表示不提供初始链，由函数内部自行初始化。
    • opts 是前面设定好的超参数和控制参数。
    • flag_stat = true 打开进度输出，flag_anim = false 关闭动画（相比原 MATLAB 代码那里默认 flag_anim = true，这是在 Octave 里为了稳定性和性能做的小改动，对推断结果没有理论影响）。910

chainer_main 内部会对每个时间点的隐状态、每个状态的均值和方差、转移矩阵等进行采样，得到一条或多条后验链 chain，用来近似整个后验分布。这个步骤在两个脚本里完全相同，只是 Octave 版本对可视化做了简化。109

5. 后验均值轨迹分析：chainer_analyze_means

  • 设置后处理参数：
    • fr = 0.25：前 25% 的样本视为 burn-in，丢弃不用。
    • dr = 2：之后每隔 2 个样本取一个，用来减弱自相关。
    • m_min = 0; m_max = 1; m_num = 25;：把状态均值的可能范围离散成 $[0,1]$ 区间上 25 个格点，用于后验平均的插值或投影。
  • 调用： [m_mod, m_red] = chainer_analyze_means(chain, fr, dr, m_min, m_max, m_num, z);
  • 输出 m_mod 可以理解为“在每一个时间点上的后验平均光子水平”，往往已被归一化到 $[0,1]$ 区间，是一个与 z 等长的一维轨迹。m_red 是某种降维或简化形式，这里没有用到。910

这一部分与 level_analysis.m 完全一致，区别在于：原脚本后面还继续做转移概率估计和漂移分析；Octave 版则转向“离散状态 + 事件检测”。109

6. 把连续水平轨迹聚类成几个离散状态

接下来是 Octave 脚本中新加的上层算法部分，用来从连续的 m_mod 中自动提取“低/中/高”强度状态：

  • 设定状态数 K = 3（可手动改成 2 或其他）。
  • 首选使用 kmeans(m_mod, K) 对 m_mod 做聚类，得到：
    • idx_raw：每个时间点所属的簇编号。
    • centers：每个簇的中心值（大致对应不同光子水平）。
  • 如果 Octave 没有统计包或 kmeans 不可用，就退化为基于分位数的简单阈值法：
    • 用三分位数 q1 = quantile(m_mod, 1/3); q2 = quantile(m_mod, 2/3); 把轨迹分为低、中、高三段。
    • 然后为每类算一个均值当作 centers。
  • 得到 centers 后，按照中心值从小到大排序，重排状态编号，使得：
    • 状态 1：最低强度（low_state）。
    • 状态 K（=3）：最高强度（high_state）。

这一段在 level_analysis.m 中没有任何对应逻辑，是 Octave 版本为 binding/bleach 检测专门设计的。其本质是：用连续的后验均值轨迹做无监督聚类，把连续空间映射为少数几个离散层级状态。910

7. 从状态序列中检测 binding 和 bleaching 事件

有了状态序列 state_seq 后，就可以做简单的“状态跳变检测”来定义事件：

  • 令 s = state_seq(:)。
  • 构造前后相邻状态：
    • s_prev = s(1:end-1);
    • s_next = s(2:end);
  • binding 事件：
    • 当某一时间点出现 s_prev == low_state 且 s_next == high_state 时，即从低强度直接跳到高强度，视为一次“分子结合”（binding）。
    • 对应的索引为 find(...) + 1，因为跳变是在第二个点发生。
  • bleach 事件：
    • 当出现 s_prev == high_state 且 s_next == low_state 时，即从高强度直接掉回低强度，视为一次“光漂白或解离”（bleaching/unbinding）。
    • 索引同样是 find(...) + 1。
  • 最后用 t(bind_idx) 和 t(bleach_idx) 把这些索引转换成真实时间点（秒），分别作为 binding_times 和 bleach_times。

这一步把 HMM 的连续后验轨迹，最终转化成了一个离散事件时间列表。在原 level_analysis.m 中，没有这类应用导向的 event detection，只停留在“描述状态分布和转移概率”的层面。109

8. 汇总所有轨迹结果并写出 CSV

  • 在 for ti = 1:n_tracks 循环内，把每条轨迹的事件信息存入 results(ti)：
    • track_id
    • binding_indices, binding_times
    • bleach_indices, bleach_times。10
  • 循环结束后，统一写出一个 CSV 文件 binding_bleach_events.csv，格式为：
    • track_index,event_type,sample_index,time_seconds
    • 其中 event_type 是 "binding" 或 "bleach"。
  • 最终在命令行打印一行提示，告诉你输出文件名。

相较之下，level_analysis.m 用的是 chainer_export 把整条 MCMC 链以 .mat 形式保存，供交互式分析和可视化使用，而不是面向“事件列表”的输出。910

总结（中文角度）

  • 从 HMM 模型和抽样的角度看： 两个脚本在 ICON HMM 的参数设定和采样过程上是同一套算法，仅在 Octave 版加入了 eps 防止除零，以及关闭动画提高稳定性，这不会改变算法本质。
  • 从整个脚本的“算法流程”看： level_analysis.m 是一个“通用 ICON HMM 演示与分析脚本”，只做单一序列的建模和一些统计性后验分析，不涉及具体事件检测逻辑。 detect_binding_bleach.m 在沿用同一 HMM 内核的基础上，增加了：多轨迹处理、对后验均值轨迹的聚类、将轨迹离散为若干强度状态、基于状态跳变的 binding/bleach 自动识别以及 CSV 输出。这是一层新的、面向应用的算法，不存在于原 MATLAB 脚本中。

如果你的目标是“在 Octave 中复用原 ICON HMM 的统计建模部分，并进一步实现自动 binding/bleaching 检测”，那么当前 detect_binding_bleach.m 的设计是符合这一思路的：核心推断算法相同，外层分析逻辑则是针对具体实验需求进行了扩展。

⁂

  1. level_analysis.m ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩

  2. detect_binding_bleach.m ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩ ↩