The problem
When we want to generate some random draws from an exotic probability distribution, the following algorithm may help:
- Generate \(u_i \sim U(0,1)\) from a uniform distribution.
- Solve \(F(x_i) = u_i\) where \(F(x)\) is the cumulative distribution function (cdf).
For most well-known distributions better methods exists, with better numerical properties especially in the tail areas. An example of such an less-known exotic distribution is the Gamma/Gompertz distribution (1).
Implementation
Here is a model I suggested to do the above:
In the above model we form a system of nonlinear equations \(G(x)=F\) in GAMS. It is a diagonal system: the Jacobian has entries on the diagonal only.
In most cases it is better to solve \(n\) single equations instead of a system of \(n\) equations. However there are a few reasons why this is not such a crazy approach as it looks:
- GAMS has a relatively high fixed cost overhead in setting up a model and calling a solver.
- A single model is a little bit more compact.
- I solved this model with CONOPT, and CONOPT will actually detect the diagonal structure and then solve the equations one by one. We can see this from some of the messages CONOPT produces.
Diagnostic messages
The log shows the following:
| C O N O P T 3 version 3.17C
Pre-triangular equations: 100 1 0 0.0000000000E+00 (After pre-processing) ** Feasible solution to a square system. |
When CONOPT receives the model it searches for triangular pieces of the model. The following graphic illustrates this:
In this picture rows are equations and columns are variables. After reordering the rows and columns we have in left-upper corner a triangular part. Those equations we can solve one-by-one. In the middle we have the remaining simultaneous equations. This is the difficult part to solve. Finally we have a post-triangular part, which we can again solve one-by-one. Note that this presolve strategy is applied for optimization models as well (in that case the Jacobian matrix is usually not a square matrix).
For our “diagonal” example all equations fall in the pre-triangular part of the model. So indeed CONOPT does not solve this as a simultaneous equation problem.
Funny message
In the listing file we found one more message:
| It may be more efficient to turn the Post-triangular equations |
Initially this reads a silly message: we have 0 post-triangular equations. However, it is the searching that is the expensive part. So even if there are eventually zero post-triangular equations, in some cases it is better to turn this search off. For our particular problem we don’t have to bother.
Results
The final results of this model look like:
| ---- 21 PARAMETER F cdf values i1 0.172, i2 0.843, i3 0.550, i4 0.301, i5 0.292, i6 0.224, i7 0.350, i8 0.856
i1 0.115, i2 0.743, i3 0.399, i4 0.205, i5 0.199, i6 0.151, i7 0.240, i8 0.767 |
Notes
As noted in the comments this particular CDF can be solved directly using:
| \[ x_i = \frac{\ln\left(\beta (1-F_i)^{-1/s}-\beta+1 \right)}{b} \] |
References
- Gamma/Gompertz distribution, https://en.wikipedia.org/wiki/Gamma/Gompertz_distribution
No comments:
Post a Comment