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

Question

在R中，在simID程序包中有一个模拟间歇性需求的函数，称为tsintermittent。

是否有等效的MATLAB？如果没有，如何在MATLAB中生成间歇性需求的随机样本？

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

Answer 1

该函数可以使用统计工具箱在MATLAB中复制。

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

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

# R Test
# 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 variance squared
% level: average of nonzero demand periods 
%
% 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