## Code prepared by Howard Hau-Chun Chuang and Rogelio Oliva ## For R version 3.0.1 http://www.r-project.org/ ## v0.01 Jan 13, 2014 ##MLE programming & OSL calculation for the Poisson-Tweedie family ##First execute the following code of the "dpoistwee" function (by Howard H. Chuang) ##The calculates the probability mass of the Poisson-Tweedie family ##The calculation is based on a double recursion algorithm by El-Shaarwi et al. (2011) ##The beginning dpoistwee=function(x,a,b,g){ probmass=rep(NA,length(x)) for(m in 1:length(x)){ p=rep(0,x[m]) if(a==0){pzero=(1-g)^b}else{pzero=exp(b*((1-g)^a-1)/a)} p[1]=b*g*pzero if(x[m]==0){prob=pzero} if(x[m]==1){prob=p[1]} if(x[m]>1){ ##Calculate the r vector beforehand r=c() r[1]=(1-a)*g for(j in 1:(x[m]-1)){ r[j+1]=((j-1+a)/(j+1))*g*r[j] } k=1 while(k<=(x[m]-1)){ j=1 sigma=0 for(j in 1:k){ sigma=sigma+p[j]*j*r[k+1-j] } p[k+1]=(b*g*p[k]+sigma)/(k+1) k=k+1 } prob=p[x[m]] } probmass[m]=log(prob) } probmass } ## ##The end ##Call the required libraries library(optimx) library(MASS) library(fitdistrplus) library(compoisson) ##Let "q" denote a vector of observed demand for a particular SKU ##We want to calculate its AIC ("PoiTw_AIC") & OSL ("PoiTw_OSL") ##The beginning ##Define the logliklihood function for the q vector PoiTwloglikf=function(par){ -sum(dpoistwee(q,a=par[1],b=par[2],g=par[3])) } ##Do a two-stage MLE ##Stage one PoiTwloglikf_1=function(par){ -sum(dpoistwee(q,a=0,b=par[1],g=par[2])) } theta0=nlminb(c(1,0.5),PoiTwloglikf_1,lower=c(0,0),upper=c(Inf,1))$par ##Stage two PoiTwfit=optim(c(0,theta0[1],theta0[2]),PoiTwloglikf) PoiTwloglik=-1*PoiTwfit$value ##Obtain the AIC value PoiTw_AIC=2*3-2*PoiTwloglik ##Now we move to OSL k=unique(q) OSLtempPoiTw=0 for(j in 1:length(k)){ qholdout=c(rep(k[j],sum(q==k[j])-1),q[q!=k[j]]) PoiTwloglikf=function(par){ -sum(dpoistwee(qholdout,a=par[1],b=par[2],g=par[3])) } ##Stage one PoiTwloglikf_1=function(par){ -sum(dpoistwee(qholdout,a=0,b=par[1],g=par[2])) } theta0=nlminb(c(1,0.5),PoiTwloglikf_1,lower=c(0,0),upper=c(Inf,1))$par ##Stage two theta=optim(c(0,theta0[1],theta0[2]),PoiTwloglikf)$par OSLtempPoiTw=OSLtempPoiTw+dpoistwee(k[j],theta[1],theta[2],theta[3])*sum(q==k[j]) } ##Obtain the OSL value PoiTw_OSL=OSLtempPoiTw ##The end ##MLE programming & OSL calculation for the CMP distribution ##First execute the following two functions "com.log.z" & "dcom.log" ##The two functions are adapted from the library "compoisson" ##The beginning com.log.z=function(lambda, nu, log.error = 0.001){ ##Perform argument checking if(lambda < 0 || nu < 0) stop("Invalid arguments, only defined for lambda >= 0, nu >= 0"); ##Initialize values z = -Inf; z.last = 0; j = 0; ## Continue until we have reached specified precision while (abs(z - z.last) > log.error & j<10000){ z.last = z; z = com.log.sum(z, j*log(lambda)-nu*com.log.factorial(j)); j = j + 1; } return (z); } dcom.log=function(x, lambda, nu, log.z = NULL){ ##Perform argument checking if(lambda < 0 || nu < 0) stop("Invalid arguments, only defined for lambda >= 0, nu >= 0"); if(x < 0 || x != floor(x)) return (0); if(is.null(log.z)) { log.z = com.log.z(lambda, nu); ##Return log pmf return((x * log(lambda) - nu * com.log.factorial(x)) - log.z); } ## ##The end ##Let "q" denote a vector of observed demand for a particular SKU ##We want to calculate its AIC ("CMP_AIC") & OSL ("CMP_OSL") ##The beginning ##Define the logliklihood function CMPloglikf=function(par){ -sum(dcom.log(q,lambda=par[1],nu=par[2])) } ##Do a two-stage MLE ##Stage one CMPloglikf_1=function(par){ -sum(dcom.log(q, lambda=mean(q),nu=par)) } nu0=nlminb(0.5,CMPloglikf_1,lower=0)$par ##Stage two CMPfit=nlminb(c(1,nu0),CMPloglikf,lower=c(1e-10,0)) CMPloglik=-1*CMPfit$objective ##Obtain the AIC value CMP_AIC=2*2-2*CMPloglik ##Now we move to OSL k=unique(q) OSLtempCMP=0 for(j in 1:length(k)){ qholdout=c(rep(k[j],sum(q==k[j])-1),q[q!=k[j]]) CMPloglikf=function(par){ -sum(dcom.log(qholdout,lambda=par[1],nu=par[2])) } ##Stage one CMPloglikf_1=function(par){ -sum(dcom.log(qholdout,lambda=mean(q),nu=par)) } nu0=nlminb(0.5,CMPloglikf_1,lower=0)$par ##Stage two CMPPar=nlminb(c(mean(qholdout),nu0),CMPloglikf,lower=c(1e-10,0))$par OSLtempCMP=OSLtempCMP+dcom.log(k[j],CMPPar[1],CMPPar[2])*sum(q==k[j]) } ##Obtain the OSL value CMP_OSL=OSLtempCMP ## ##The end ##MLE programming & OSL calculation for the NB distribution ##Let "q" denote a vector of observed demand for a particular SKU ##We want to calculate its AIC ("NB_AIC") & OSL ("NB_OSL") ##The beginning NBfit=fitdistr(q,"Negative Binomial") ##Obtain the AIC value NB_AIC=2*2-2*NBfit$loglik ##Now we move to OSL k=unique(q) OSLtempNB=0 for(j in 1:length(k)){ qholdout=c(rep(k[j],sum(q==k[j])-1),q[q!=k[j]]) NBPar=fitdistr(qholdout,"Negative Binomial")$estimate OSLtempNB=OSLtempNB+dnbinom(k[j],size=NBPar[1],mu=NBPar[2],log=T)*sum(q==k[j]) } ##Obtain the OSL value NB_OSL=OSLtempNB ## ##The end