Docs

docs/src/customize.md

@@ -183,7 +183,7 @@ penalize(d::GammaExpert, p) = (p[1]-1)*log(d.k) - d.k/p[2] + (p[3]-1)*log(d.θ)
 
 In the function `penalize`, the hyperparameters are given as a vector `p`.
 We assume the shape parameter `k` of the Gamma expert follows a Gamma prior distribution with shape `p[1]` and scale `p[2]`. Analogously, the scale parameter
-`θ` of the Gamma expert follows a Gamma prior distribution with shape `p[1]` and scale `p[2]` of the Gamma expert follows a Gamma prior distribution with shape `p[3]` and scale `p[3]`. The `penalize` function calculates the prior logpdf of the
+`θ` of the Gamma expert follows a Gamma prior distribution with shape `p[3]` and scale `p[4]`. The `penalize` function calculates the prior logpdf of the
 parameters, excluding the constant terms since they are irrelevant to the EM algorithm.
 
 The functions `penalty_init` initializes some default hyperparameters, if they
@@ -192,17 +192,110 @@ which poses no penalty, as plugging ion `p = [1.0 Inf 1.0 Inf]` into the
 `penalize` function yields zero. Still, it is recommended to always use some penalty
 in application to avoid the issue of spurious models mentioned above.
 
+Note that the penalty term should be subtracted from the loglikelihood
+when implementing the M-step of the EM algorithm (see below).
 
+**Miscellaneous functions**: There are a number of miscellaneous functions implemented for the expert function as shown below. They are mainly used in 
+predictive functions such as `predict_mean_prior`. It is recommended to code them
+as well when adding a new expert function, but (some or all of) these functions can be optional if
+the user only wants to fit an LRMoE model.
 
+```julia
+## statistics
+mean(d::GammaExpert) = mean(Distributions.Gamma(d.k, d.θ))
+var(d::GammaExpert) = var(Distributions.Gamma(d.k, d.θ))
+quantile(d::GammaExpert, p) = quantile(Distributions.Gamma(d.k, d.θ), p)
+lev(d::GammaExpert, u) = d.θ*d.k*gamma_inc(float(d.k+1), u/d.θ, 0)[1] + u*(1-gamma_inc(float(d.k), u/d.θ,0)[1])
+excess(d::GammaExpert, u) = mean(d) - lev(d, u)
+```
 
 
 
 
 
+### M-Step: exact observation
 
-### Maximizing loglikelihood: exact observation
-
-The next step is to derive the loglikelihood of a Gamma distribution, which is given in Section 3.2 of [Fung et al. (2019)](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3740061). Notice that the cited paper considers the issue of data truncation and censoring,
+The next step is to implement the M-step of the EM algorithm, which is given in Section 3.2 of [Fung et al. (2019)](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3740061). Notice that the cited paper considers the issue of data truncation and censoring,
 but as a first step, we will just assume all data are observed exactly. Data truncation and censoring are dealt with
 afterwards.
 
+The M-step with exact observation is carried out by the function
+`EM_M_expert_exact` in the source code. Note that the function name and arguments
+should not be altered, since it is referenced in the main fitting function of the
+package.
+
+In short, the goal is to maximize the objective function  `z_e_obs` * LL, where LL is
+the loglikelihood of the expert `d` when observing a vector `ye`. If the function argument
+`penalty` is true, then a penalty term given by hyperparameters `pen_params_jk` (see also above) is subtracted from the objective function. The function argument
+`expert_ll_pos` is a legacy and irrelevant to the M-step when observations `ye` are exact.
+
+For Gamma expert, we are maximizing with respect to ``k`` and ``\lambda`` (without penalty) the objective function `z_e_obs` multiplied by
+
+``LL = \sum_{i = 1}^{n} (k-1)\log(y_i) - y_i/k - \log(\Gamma(k)) - k\log(\lambda)``.
+
+The optimization is done by calling `minimizer` in `Optim.jl`, which is ommited in this document. Penalty can also be added in the objective function if
+the function argument `penalty` is set to `true`.
+
+Finally, the `EM_M_expert_exact` returns a Gamma expert object containing
+the updated parameters, thus completing the M-step.
+
+### M-Step: censored and truncated observation
+
+One important feature of this package is to deal with censored and truncated data,
+which is common both in and outside actuarial contexts. In this case, the
+loglikelihood is incomplete, in the sense that ``y_i`` above is not known exactly,
+leading to an additional E-step before the M-step.
+
+In particular, the M-step with data censoring and truncation is done by the 
+function `EM_M_step`, which takes a few more arguments compared with 
+`EM_M_expert_exact`. Again, the function name and arguments should not be altered.
+
+Recall that we assume the obsercations are truncated between `tl` and `tu`, as well as censored between `yl` and `yu`. These lower and upper bounds are needed for the
+M-step.
+
+Censored data points are in the observed dataset. For all the censored observation `y`, we don't know their exact values, but only know that they are between 
+`yl` and `yu`. Hence, for those observations, we should take the expected value
+of the loglikelihood, further normalizing them by the probability of falling within
+this interval (which is equal to `exp.(-expert_ll_pos)`).
+
+For the Gamma expert, we refer to the equation of ``LL`` above: the uncertain terms
+are ``\log(y_i)`` and ``y_i``, so we should compute their conditional expectation,
+given that the exact observation is between `yl` and `yu`. This corresponds to
+some numerical integration procedures in the source code. Note that the conditional
+expectation should be computed with respect to the probability measure
+implied by the old (i.e. not yet updated) parameters.
+
+A similar remark can be made about truncated observations, which are not present
+in the dataset due to truncation. Those are either smaller than the lower bound
+of truncation `tl`, or larger than the upper bound of truncation `tu`. In addition,
+the function argument `k_e` is equal to the expected number of lost observations
+due to truncation, which should also be multiplied to the conditional expectation
+of loglikelihood.
+
+To be more specific, we also numerically integrate ``\log(y_i)`` and ``y_i`` in ``LL``, but from `0` to `tl` and from `tu` to `Inf`. Then, we normalize it by the probability `exp.(-expert_tn_bar_pos)` to obtain the conditional expectation.
+Finally, the loglikelihood is further multiplied by `k_e` to account for all
+unobserved data points due to truncation.
+
+Finally, the loglikelihood of truncated observations is added to the objective function for optimization.
+
+## Notes and tips
+
+There are quite a few things to note when adding a customized expert function.
+
+* Exact observations: If the problem at hand only concerns exact observations,
+  then the user can write the `EM_M_exact` function only.
+
+* Speeding up numerical integration: In `EM_M_step`, numerical integration
+  for certain expert functions can take a long time. In the `gamma.jl` example,
+  the author has first looked up for unique lower and upper bounds, and only integrate using those pairs. Then, the numerical integration is mapped to the original dataset. Users are advised to also take this approach, which can be done by conveniently copying and modifying the source code of `gamma.jl`.
+
+* Optimization constraits: Some expert functions pose constrants on the parameters,
+  e.g. the parameters of Gamma experts should be positive. However, constrained optimization can be computationally slow. In `gamma.jl`, we have used a log transformation, so that the optimization procedure `Optim.optimize` is unconstrained. Also, specifying a small search interval speeds things up.
+
+* Zero inflation: For zero-inflated expert functions, the M-step should only
+  take in positive observations for continuous experts, but all non-negative observations for discrete experts. See `gamma.jl` and `poisson.jl` for a comparison.
+
+For further questions and tips when customizing expert functions, please contact the author on github.
+
+
+