我已经实现了类似的东西来计算气流中粒子的速率/浓度。它们来自随机事件(可能是泊松分布),我想知道每次粒子的平均速率。我的方法如下(取自我的代码的文档字符串):
事件时间戳被放置在一个固定的最大大小的缓冲区中。
每当请求浓度估计时,所有时间戳
过滤掉设置的最大年龄。如果剩余数量
据报道,最近的时间戳低于某个阈值
浓度为零。否则,计算浓度
剩余事件数除以时间差
在缓冲区中最早的时间戳和现在之间。
我将附上下面的 Python 实现以供参考。它是一个更大项目的一部分,所以我不得不稍微修改它以摆脱对外部代码的引用:
particle_concentration_estimator_params.py
#!/usr/bin/env python3
"""This module implements a definition of parameters for the particle
concentration estimator.
"""
__author__ = "bup"
__email__ = "bup@swisens.ch"
__copyright__ = "Copyright 2021, Swisens AG"
__license__ = "GNU General Public License Version 3 (GPLv3)"
__all__ = ['ParticleConcentrationEstimatorParams']
import dataclasses
@dataclasses.dataclass
class ParticleConcentrationEstimatorParams:
"""This provides storage for the parameters used for the particle
concentration estimator.
"""
timestamp_buffer_max_size: int
"""Maximum size of the buffer that is used to keep track of event
timestamps. The size of this buffer mainly affects the filtering
of the reported data.
Unit: - (count)
"""
timestamp_buffer_max_age: float
"""Maximum age of events in the timestamp buffer which are
considered for the concentration calculation. This value is a
tradeoff between a smooth filtered value and the dynamic response
to a changed concentration.
Unit: s
"""
min_number_of_timestamps: int
"""Minimum number of timestamps to use for the concentration
estimation. If less timestamps are available, the concentration is
reported as zero.
Unit: - (count)
"""
particle_concentration_estimator.py
#!/usr/bin/env python3
"""This module implements the particle concentration estimation.
"""
__author__ = "bup"
__email__ = "bup@swisens.ch"
__copyright__ = "Copyright 2021, Swisens AG"
__license__ = "GNU General Public License Version 3 (GPLv3)"
__all__ = ['ParticleConcentrationEstimator']
import logging
import time
from typing import Optional
import numpy as np
from .particle_concentration_estimator_params import ParticleConcentrationEstimatorParams
logger = logging.getLogger(__name__)
class ParticleConcentrationEstimator:
"""An object of this class implements the Poleno particle
concentration estimator. Particle concentration is basically just
a number of particles per time unit. But since the particle events
arrive irregularly, there are various ways to filter the result, to
avoid too much noise especially when the concentration is low. This
class implements the following approach:
Event timestamps are placed in a buffer of a fixed maximum size.
Whenever a concentration estimate is requested, all timestamps up to
a set maximum age are filtered out. If the remaining number of
recent timestamps is below some threshold, the reported
concentration is zero. Otherwise, the concentration is calculated
as the number of remaining events divided by the time difference
between the oldest timestamp in the buffer and now.
"""
def __init__(self, params: ParticleConcentrationEstimatorParams):
"""Initializes the object with no events.
Args:
est_params: Initialized PolenoParams object which includes
information describing the estimator's behaviour.
"""
self.params = params
"""Stored params for the object."""
n_rows = self.params.timestamp_buffer_max_size
self._rb = np.full((n_rows, 2), -1e12)
self._rb_wp = 0
self._num_timestamps = 0
self._concentration_value = 0.0
self._concentration_value_no_mult = 0.0
def tick(self, now: float) -> float:
"""Recalculates the current concentration value.
Args:
now: Current timestamp to use to filter out old entries
in the buffer.
Returns:
The updated concentration value, which is also returned
using the concentration attribute.
"""
min_ts = now - self.params.timestamp_buffer_max_age
min_num = self.params.min_number_of_timestamps
used_rows = self._rb[:, 0] >= min_ts
filt_ts = self._rb[used_rows]
num_ts = round(np.sum(filt_ts[:, 1]))
self._num_timestamps = num_ts
num_ts_no_mult = round(np.sum(used_rows))
if num_ts < min_num:
self._concentration_value = 0.0
self._concentration_value_no_mult = 0.0
else:
t_diff = now - np.min(filt_ts[:, 0])
if t_diff >= 1e-3:
# Do not change the reported value if all events in the
# buffer have the same timestamp.
self._concentration_value = num_ts / t_diff
self._concentration_value_no_mult = num_ts_no_mult / t_diff
return self._concentration_value
def got_timestamp(self,
ts: Optional[float] = None,
multiplier: float = 1.0) -> None:
"""Passes in the most recent event timestamp. Timestamps need
not be ordered.
Calling this method does not immediately update the
concentration value, this is deferred to the tick() method.
Args:
ts: Event timestamp to use. If None, the current time is
used.
multiplier: Optional multiplier, which makes ts to be
counted that many times.
"""
if ts is None:
ts = time.time()
self._rb[self._rb_wp] = (ts, float(multiplier))
self._rb_wp = (self._rb_wp + 1) % self._rb.shape[0]
self._num_timestamps += round(multiplier)
@property
def concentration(self) -> float:
"""The calculated concentration value.
Unit: 1/s
"""
return self._concentration_value
@property
def concentration_no_mult(self) -> float:
"""The calculated concentration value without taking into
account the timestamp multipliers, i.e. as if all timestamps
were given with the default multiplier of 1.
Unit: 1/s
"""
return self._concentration_value_no_mult
@property
def num_timestamps(self) -> int:
"""Returns the number of timestamps which currently are in
the internal buffer.
"""
return self._num_timestamps