如何在MATLAB中模拟间歇性需求？

Question

背景：间歇性需求是指（需求）时间序列数据，其特征在于存在需求为零的期间。这些非零需求通常很小但变化很大。间歇性需求有许多应用，并给许多传统的预测和控制系统带来挑战。

如何在MATLAB中生成间歇性需求的随机样本？

在R的tsintermittent程序包中，有一个用于模拟间歇性需求的函数simID。
我不知道有任何内置的MATLAB函数可以完成此操作（我非常怀疑它们是否存在）。

相关参数包括可变性和间歇性。

Answer 1

可以使用nbinrnd和使用Statistics工具箱在MATLAB中复制该函数binornd。

语法[d] = simID(n, obs, adi, cv2, level)与

n：创建
obs的需求流数量：
adi的模拟时间段数：平均需求间隔（时间段）
cv2：变异系数平方
水平：非零需求的平均值
d：'n'x'obs'矩阵，具有指定的需求流

这明确地假设非零需求到达遵循伯努利分布，并且非零需求达到负二项式分布。^[1]

_{[1] Petropoulos F.，Makridakis S.，Assimakopoulos V.和Nikolopoulos K.（2014）“需求预测中的“课程的马匹””，《欧洲运筹学杂志》，第1卷。237，第1号，第152-163页。}

% Test in MATLAB 
D = simID(1,50000,2.5,2.5,3);

mean(D)                          % 1.2089
std(D)                           % 3.3913

# Test in R
# install.packages('tsintermittent')
library('tsintermittent')

D = simID(1,50000,2.5,2.5,3)
mean(D$ts.1)                      # 1.19632
sd(D$ts.1)                        # 3.261225

由于没有相同的随机种子，所以值不完全一致。

MATLAB代码：（已通过R2019a测试）

function [d] = simID(n, obs, adi, cv2, level)
% n: number of demand streams to create
% obs: number of time periods to simulate
% adi: average demand interval
% cv2: coefficient of variation squared
% level: average of nonzero demand 
%
% d: 'n' x 'obs' matrix with the specified demand streams
%
% If 'level' empty, defaults to sample from Uniform(10,100). 
% 
% MATLAB implementation of simID() from 'tsintermittent' package for R
% https://www.rdocumentation.org/packages/tsintermittent/versions/1.9/topics/simID
% https://CRAN.R-project.org/package=tsintermittent
%
% This simulator assumes that non-zero demand arrivals follow a bernoulli distribution 
% and the non-zero demands a negative binomial distribution. 
% Petropoulos F., Makridakis S., Assimakopoulos V. & Nikolopoulos K. (2014) 
% "'Horses for Courses' in demand forecasting", 
% European Journal of Operational Research, Vol. 237, No. 1, pp. 152-163

% Input Error Checking ****************************************************
narginchk(3,5)
if isempty(level), level = 10 + 90*rand(); end
if cv2<=0, error('cv2 must be positive'), end
if n<0 || obs<0 || adi < 0, error('n, obs, and adi must be positive'), end
% End (Input Error Checking) **********************************************

m = level - 1;
if ( cv2 ~= 0)
    p = (m / ( cv2 * ( (m+1)^2) ) ); % calculates the p for negative binomial function
    r = m * p / (1 - p); % calculates the r for the negative binomial funcion
    d = binornd(1, 1 / adi, n, obs) .* (nbinrnd(r, p, n, obs) + 1);
else
    d = binornd(1, 1/adi, n, obs) .* repmat(m + 1, n, obs);
end

end