|
ContactCenters V. 0.9.9. |
||||||||
PREV CLASS NEXT CLASS | FRAMES NO FRAMES | ||||||||
SUMMARY: NESTED | FIELD | CONSTR | METHOD | DETAIL: FIELD | CONSTR | METHOD |
java.lang.Objectumontreal.iro.lecuyer.contactcenters.contact.GammaParameterEstimator
public class GammaParameterEstimator
This class implements the parameter estimation for the doubly Gamma-Poisson process.
The rate Ti, j of this process consists of three multiplicative contributions: (i)
the deterministic piece-wise constant rate λi, (ii) the busyness factor for
the day, βj, (iii) the business factor for the sub-period of the day, Bi, j.
The input data are counts observed for I sub-periods of the day during J days.
We assume that the rate follows
Ti, j = λiβjBi, j,∀i = 1…I, j = 1…J
and input data
Yi, j∼Poisson(Ti, j) follow the Poisson distribution
conditional on the rates. The busyness factor for the day follows Gamma distribution with
parameters Q and Q. The busyness factor follows the Gamma distribution with parameters
R and R. The busyness factors are assumed to be independent across
days and sub-intervals of the day. The class implements two estimators of the process parameters:
the Moment Matching Estimator (MME) in the method getMMEdoublyGamma()
and the Maximum
Likelihood Estimator (MLE) in the method getMLEdoublyGamma()
. The MME is a
suboptimal, but simple and fast estimator based on matching the theoretical means, variances
and covariances of the process with their empirical counterparts. The MLE is a statistically
optimal estimator with greater accuracy than the MME, but it is less computationally efficient.
The MLE is implemented using the stochastic trust-region Gauss-Newton algorithm. The MLE estimator first calls the MME estimator
and uses it as a starting point for the optimization. The startGradUncTrustRegion(double[])
method uses
common random numbers to track changes in the cost function and call of this method with option
"CostOnly" must, in general, be always preceded by the call of this method without this option
in order to provide meaningful results.
Field Summary | |
---|---|
double[][] |
x
|
Constructor Summary | |
---|---|
GammaParameterEstimator(int[][] data,
int N,
int P)
Constructs a new estimator object with a given set of input data |
|
GammaParameterEstimator(int[][] data,
int N,
int P,
int M)
Constructs a new estimator object with a given set of input data and default seed for the Gamma random variable generators. |
|
GammaParameterEstimator(int[][] data,
int N,
int P,
int M,
long[] Seed)
Constructs a new estimator object with a given set of input data and a user defined seed for the Gamma random variable generators. |
Method Summary | |
---|---|
void |
estimateNortaRateParamsStochasticRootFinding()
Estimates the parameters of the Gamma-Poisson NORTA model for rates using stochastic root finding approach. |
double[] |
getLikelihoodDerivativesDoublyGamma(double[] Lam,
double R,
double Q,
String OutType)
Calculates the values of the log-likelihood function and its derivatives for the doubly Gamma-Poisson arrival process model. |
double[][] |
getLikelihoodDerivativesDoublyGammaSpline(double[] Lam,
double[] R,
double Q,
String OutType)
Calculates the values of the log-likelihood function and its derivatives for the doubly Gamma-Poisson arrival process model. |
double |
getLogNegBinDer(int[] X,
double Xmean,
double r)
Calculates the derivative of the log-likelihood function for the Negative binomial distribution. |
double[] |
getMLEdoublyGamma()
Estimates the parameters of a doubly Gamma Poisson-Gamma arrival process that has both busyness factor for the day and the busyness factor for the sub-period of the day, both following the Gamma distribution, from the number of arrivals in the array arrivals using maximum likelihood approach. |
double[] |
getMLEdoublyGammaSpline()
Estimates the parameters of a doubly Gamma Poisson-Gamma arrival process that has both busyness factor for the day and the busyness factor for the sub-period of the day, both following the Gamma distribution, from the number of arrivals in the array arrivals using maximum likelihood approach. |
double[] |
getMMEdoublyGamma()
Estimates the parameters of a doubly Gamma Poisson-Gamma arrival process that has both busyness factor for the day and the busyness factor for the sub-period of the day, both following the Gamma distribution, from the number of arrivals in the array arrivals using method of moments. |
double[] |
getMMEdoublyGammaGeneral()
Estimates the parameters of a doubly Gamma Poisson-Gamma arrival process that has both busyness factor for the day and the busyness factor for the sub-period of the day, both following the Gamma distribution, from the number of arrivals in the array arrivals using method of moments. |
double[] |
getNegBinMLE(int[] X,
double tole)
Calculates the MLEs of parameters of the negative binomial distribution. |
double[] |
getNortaRateGammaScale()
Returns the estimated vector of λG parameters of the Gamma distribution in the compound Gamma-Poisson NORTA model for rates. |
double[] |
getNortaRateGammaShape()
Returns the estimated vector of αG parameters of the Gamma distribution in the compound Gamma-Poisson NORTA model for rates. |
double[][] |
getNortaRateGaussCorr()
Returns the estimated copula correlation matrix for the Gamma-Poisson NORTA model for rates. |
double[][] |
getNortaRateGaussCorrCorrected()
Returns the estimated and corrected copula correlation matrix for the Gamma-Poisson NORTA model for rates. |
double[] |
getNortaRateGaussCorrFitGeneralLinear()
Fits the general linear model rj = abj + c to the estimated copula correlation matrix for the Gamma-Poisson NORTA model for rates. |
double |
getNortaRateGaussCorrFitMarkovSingleRho()
Fits the single ρ Markov linear model rj = bj to the estimated copula correlation matrix for the Gamma-Poisson NORTA model for rates. |
double |
getNortaRhoStochasticRootFinding(double rhoTarget,
double[][] NegBinParams,
double rhoInit)
Solves the problem of fitting the NORTA correlation coefficient to the empirical Spearman correlation coefficient of counts in the Gamma-Poisson copula model using stochastic root finding approach. |
double[][] |
getOptimizationTrace()
Returns the the matrix of the optimization trace containing the evolution of parameters during optimization iterations. |
double[] |
getPolyakAverage()
Calculates the estimates of the parameters based on the stochastic optimization framework described by B T Polyak and A B Juditsky in "Acceleration of Stochastic Approximation by Averaging", SIAM J Control Optim. |
void |
initNortaStochRootFinding()
|
void |
initNortaStochRootFinding(long[] Seed)
|
void |
setMaxIter(int maxit)
|
void |
setMovingWindowSize(int movWindSize)
Sets the number of sub-periods over which to average the MME estimate in getMMEdoublyGammaGeneral. |
void |
setSmoothingLambda(double lambda)
Sets the smoothing parameter for the use with the smoothing spline doubly gamma penalized MLE. |
void |
setTrustRegionAnnealingPwr(double Value)
Sets the power law in the step size stochastic approximation annealing sequence. |
void |
setTrustRegionInitBoundary(double Value)
Sets the initial size of the trust region (initial value of the step size regulator). |
void |
setTrustRegionMaxIterations(int Value)
Sets the maximum number of iterations in the trust region optimization algorithm |
void |
setTrustRegionMultipliers(double IncreaseMultiplier,
double DecreaseMultiplier)
Sets the multipliers that control the rate at which the step size regulator is increased ( g1 )
or decreased (g0 ). |
void |
setTrustRegionNoiseAttenuator(double Value)
Sets the value of the noise attenuator for the trust region approach. |
void |
setTrustRegionQualBounds(double LowQualBound,
double HighQualBound)
Sets the lower and higher bounds on the quality metrics of the trust region. |
void |
setTrustRegionTol(double Value)
Sets the threshold for the tolerance of the trust region optimization algorithm. |
double[] |
startGradUncTrustRegionSpline(double[] SolInit)
Implements the trust region Gauss-Newton maximizer of the likelihood of the doubly Gamma-Poisson process. |
Methods inherited from class java.lang.Object |
---|
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait |
Field Detail |
---|
public double[][] x
Constructor Detail |
---|
public GammaParameterEstimator(int[][] data, int N, int P)
data
- the matrix of input data Yi, j with N rows corresponding to N
observations and P columns corresponding to P sub-periods in the day.N
- the number of observationsP
- the number of sub-periods in the daypublic GammaParameterEstimator(int[][] data, int N, int P, int M)
data
- the matrix of input data Yi, j with N rows corresponding to N
observations and P columns corresponding to P sub-periods in the day.N
- the number of observationsP
- the number of sub-periods in the dayM
- the number of Monte-Carlo samples used in the evaluation of
stochastic derivatives and the cost functionpublic GammaParameterEstimator(int[][] data, int N, int P, int M, long[] Seed)
data
- the matrix of input data Yi, j with N rows corresponding to N
observations and P columns corresponding to P sub-periods in the day.N
- the number of observationsP
- the number of sub-periods in the dayM
- the number of Monte-Carlo samples used in the evaluation of
stochastic derivatives and the cost functionSeed
- is the vector of 6 integers. The first 3 values of the seed
must all be less than
m1 = 4294967087, and not all 0; and the last 3
values must all be less than
m2 = 4294944443, and not all 0.Method Detail |
---|
public void initNortaStochRootFinding(long[] Seed)
public void initNortaStochRootFinding()
public double[] getMMEdoublyGamma()
public double[] getMMEdoublyGammaGeneral()
public double[] getMLEdoublyGamma()
startGradUncTrustRegion(double[])
and the stochastic approximation of the likelihood function and its derivatives implemented
in getLikelihoodDerivativesDoublyGamma(double[], double, double, java.lang.String)
in order to obtain the MLE. Before
launching the Gauss-Newton optimizer, this algorithm calls the getMMEdoublyGamma()
in order
to get good initialization for the values of parameters. The output is formatted the
same way as it is done for the getMMEdoublyGamma()
.
If we follow the notation introduced for PoissonGammaArrivalProcess, then this method
estimates
αG, p,
λG, p and the daily gamma
busyness parameter. It is assumed that
αG, p is a vector of
distinct values while all the entries of
λG, p are the same and equal to R.
Thus the returned array of 2P + 1 elements contains
(αG, 0,…, αG, P-1, λG, 0,…, λG, P-1, β0.
public double[] getMLEdoublyGammaSpline()
startGradUncTrustRegion(double[])
and the stochastic approximation of the likelihood function and its derivatives implemented
in getLikelihoodDerivativesDoublyGamma(double[], double, double, java.lang.String)
in order to obtain the MLE. Before
launching the Gauss-Newton optimizer, this algorithm calls the getMMEdoublyGamma()
in order
to get good initialization for the values of parameters. The output is formatted the
same way as it is done for the getMMEdoublyGamma()
.
If we follow the notation introduced for PoissonGammaArrivalProcess, then this method
estimates
αG, p,
λG, p and the daily gamma
busyness parameter. It is assumed that
αG, p is a vector of
distinct values while all the entries of
λG, p are the same and equal to R.
Thus the returned array of 2P + 1 elements contains
(αG, 0,…, αG, P-1, λG, 0,…, λG, P-1, β0.
public double[] getLikelihoodDerivativesDoublyGamma(double[] Lam, double R, double Q, String OutType)
For the doubly Gamma-Poisson arrival process model the log-likelihood function
does not admit any closed form as it contains
an integral over βj, the daily busyness factor. The integral cannot be treated
analytically. This integral is thus treated numerically via the Monte-Carlo
approach. The Monte-Calro approach uses numSamples
samples for the
evaluation of the integral. The Gamma distribution with both parameters equal to
Q is taken as the proposal distribution. The cost function and all the derivatives
are thus approximated stochastically. There is an option to return only the
value of the cost function (the log-likelihood function) by assigning to OutType
the string "CostOnly". In this case the algorithm uses the random numbers common with
the ones that were used to evaluate the derivatives last time the function was called.
The importance sampling is used to compensate for possible change of integration measure
(that may arise due to change in the value of Q between calls). This feature is used
by the trust region algorithm for determining the quality of the trust region and for
regulating the step size of the optimization. The common random numbers approach reduces
the effects of noise on the step size regulation.
The returned array of 2*(P + 2) + 1 elements contains
1) The value of the cost function 2) P values of first order derivatives for deterministic rates 3) Derivative with respect to $R$, sub-period Gamma rate 4) Derivative with respect to $Q$, daily Gamma rate 5) P values of second order derivatives for deterministic rates 6) Second order derivative with respect to $R$, sub-period Gamma rate 7) Second order derivative with respect to $Q$, daily Gamma rate
Lam
- initial values of base rates.R
- initial value of the Gamma distribution parameter for the sub-period of the day business factor.Q
- initial value of the Gamma distribution parameter for the daily business factor.OutType
- the string controlling output options
public double[][] getLikelihoodDerivativesDoublyGammaSpline(double[] Lam, double[] R, double Q, String OutType)
For the doubly Gamma-Poisson arrival process model the log-likelihood function
does not admit any closed form as it contains
an integral over βj, the daily busyness factor. The integral cannot be treated
analytically. This integral is thus treated numerically via the Monte-Carlo
approach. The Monte-Calro approach uses numSamples
samples for the
evaluation of the integral. The Gamma distribution with both parameters equal to
Q is taken as the proposal distribution. The cost function and all the derivatives
are thus approximated stochastically. There is an option to return only the
value of the cost function (the log-likelihood function) by assigning to OutType
the string "CostOnly". In this case the algorithm uses the random numbers common with
the ones that were used to evaluate the derivatives last time the function was called.
The importance sampling is used to compensate for possible change of integration measure
(that may arise due to change in the value of Q between calls). This feature is used
by the trust region algorithm for determining the quality of the trust region and for
regulating the step size of the optimization. The common random numbers approach reduces
the effects of noise on the step size regulation.
The returned 2D array Output of P×7 elements contains
1) Output[0][0] The value of the cost function 2) Output[:][1] (second column) P values of first order derivatives with respect to deterministic rates 3) Output[:][2] P values of second order derivatives with respect to deterministic rates 4) Output[:][3] P values of first order derivatives with respect to $R_p$, sub-period Gamma rate 5) Output[:][4] P values of second order derivatives with respect to $R_p$, sub-period Gamma rate 6) Output[0][5] Derivative with respect to $Q$, daily Gamma rate 7) Output[0][6] Second order derivative with respect to $Q$, daily Gamma rate
Lam
- initial values of base rates.R
- initial value of the Gamma distribution parameter for the sub-period of the day business factor.Q
- initial value of the Gamma distribution parameter for the daily business factor.OutType
- the string controlling output options
public double[] startGradUncTrustRegionSpline(double[] SolInit)
eta
, which is a positive
constant smaller than one.
In the case eta
is smaller than 1, only a portion of the gradient is
used during each iteration. As the number of Monte-Carlo samples used to approximate the
derivatives reduces, the (absolute) value of the noise attenuator must also be reduced.
The second option is to use the stochastic approximation approach with the step size
reduction sequence of the form
t-α, where t is
the iteration number and α is the number between 1/2 and 1. In our implementation
we have parameter pwr
with default value 5/6 that determines the speed of
decay of the step sequence and is equivalent to α.
SolInit
- the value of the solution at first iteration
public void setTrustRegionMaxIterations(int Value)
Value
- maximum number of iterations in the trust region optimization algorithmpublic void setTrustRegionQualBounds(double LowQualBound, double HighQualBound)
c0
then the size of the trust region
(the step size regulator) will be reduced. If the quality metric has value higher than
c1
then the size of the trust region (the step size regulator) will be increased.
LowQualBound
- is the lower boundHighQualBound
- is the higher boundpublic void setTrustRegionMultipliers(double IncreaseMultiplier, double DecreaseMultiplier)
g1
)
or decreased (g0
). If the quality metric has value lower than c0
then the size of
the trust region (the step size regulator) is multiplied by g0
. If the quality
metric has value higher than c1
then the size of the trust region
(the step size regulator) will be multiplied by g1
. IncreaseMultiplier must be greater than 1,
DecreaseMultiplier must be between 0 and 1.
IncreaseMultiplier
- DecreaseMultiplier
- public void setTrustRegionTol(double Value)
Value
- the value of the stopping criterionpublic void setTrustRegionNoiseAttenuator(double Value)
Value
- the value of the noise attenuatorpublic void setTrustRegionInitBoundary(double Value)
Value
- the initial value of the step size regulator.public void setTrustRegionAnnealingPwr(double Value)
Value
- public void setSmoothingLambda(double lambda)
lambda
- smoothing parameter in the interval [0, 1]. When this parameter is equal to 1, the
smoothing spline is not used and shape parameters are assumed to be independent
across sub-periods. Default value 0.95.public void setMovingWindowSize(int movWindSize)
movWindSize
- number of sub-periods over which to average the MME estimate.public void setMaxIter(int maxit)
public double[][] getOptimizationTrace()
public double[] getPolyakAverage()
Rinit
0.1, eta
0.3 and 5/6
.
public void estimateNortaRateParamsStochasticRootFinding()
yGaussCorr
, the estimated base rates in Qout
and the parameter of the Gamma distribution in LamOut
.
These quantities can be accessed using methods getNortaRateGaussCorr()
,
getNortaRateGammaShape()
and getNortaRateGammaScale()
.
The parameters of the Gamma distribution are estimated using method getNegBinMLE(int[], double)
.
The entries of the copula correlation matrix are estimated using
method getNortaRhoStochasticRootFinding(double, double[][], double)
.
public double[] getNegBinMLE(int[] X, double tole)
X
- vector of observed countstole
- parameter search tolerance for the binary search
public double getLogNegBinDer(int[] X, double Xmean, double r)
X
- data valuesXmean
- mean of the data valuesr
- parameter of the negative binomial distribution (number of failures)
public double getNortaRhoStochasticRootFinding(double rhoTarget, double[][] NegBinParams, double rhoInit)
rhoTarget
- the empirically observed Spearman correlation coefficient
of counts in the Gamma-Poisson copula modelNegBinParams
- estimated parameters of the marginal distribution,
which is Negative Binomial in the case of Gamma-Poisson copula modelrhoInit
- initial value of the NORTA correlation coefficient. Typically,
rhoInit = rhoTarget
public double[] getNortaRateGammaShape()
PiecewiseConstantPoissonArrivalProcess
so that the base rate
in subperiod p is equal to
αG, p/λG, p.
public double[] getNortaRateGammaScale()
PiecewiseConstantPoissonArrivalProcess
so that the base rate
in subperiod p is equal to
αG, p/λG, p.
public double[][] getNortaRateGaussCorr()
public double[][] getNortaRateGaussCorrCorrected()
public double getNortaRateGaussCorrFitMarkovSingleRho()
public double[] getNortaRateGaussCorrFitGeneralLinear()
|
ContactCenters V. 0.9.9. |
||||||||
PREV CLASS NEXT CLASS | FRAMES NO FRAMES | ||||||||
SUMMARY: NESTED | FIELD | CONSTR | METHOD | DETAIL: FIELD | CONSTR | METHOD |