Model Selection Guided Tutorials

One way that you can provide a smooth estimate of a density without reference to a parametric family is through a so-called kernel density estimate. A kernel density estimate of a probability density function (f(x)) has the following form:
$$ f_n(x) = frac{1}{nb} sum_{i=1}^n kleft(frac{x-X_i}{b}right).$$ In this expression, (X_1, ldots, X_n) is our random sample of (n) observations, the positive constant (b) is known as the bandwidth and the function ( k(cdot ) ) is called a kernel function. The following are standard choices of the kernel function:

• uniform kernel, ( k(y) = frac{1}{2} I(|y| le 1) )
• triangular kernel, (k(y) = (1-|y|)times I(|y| le 1))
• Epanechnikov kernel, (k(y) = frac{3}{4}(1-y^2) times I(|y| le 1))
• Gaussian kernel (k(y) = phi(y)), where (phi(cdot)) is the standard normal density function.

For some situations, you can also use kernel methods to give a smooth approximation of the distribution function as follows:
$$ hat{F}_n(x) = frac{1}{n} sum_{i=1}^n Kleft(frac{x-X_i}{b}right).$$ Here, ( hat{F}_n(x)) is known as the kernel density estimator of a distribution function. The function (K ) is a probability distribution function associated with the kernel density (k). To illustrate, for the uniform kernel, we have
$$ K(y) = begin{cases}
0 & y<-1\
frac{y+1}{2}& -1 le y < 1 \
1 & y ge 1 .\ end{cases} $$ To solve questions using kernel smoothing, you can do the following:

1. From the problem, identify the data points (X_i), the kernel function to be applied, and the bandwidth (b) .
2. Rescale the observations about a point (x). Call the rescaled observations ( y_i = (x – X_i)/b ).
3. Apply the kernel density ( ( k(y_i) ) ) or distribution function ( ( K(y_i) ) ) , as appropriate. Kernels such as uniform, triangular and Epanechinikov have a limited domain, so make sure that you account for this feature when you evaluate the kernel function.
4. Then, take the average. For the density estimator, divide this average by (b). (Not needed for the distribution function.)

Suppose that you have right-censored data and you wish to estimate the distribution or survival function without reference to a parametric model. Two options are available, one estimator due to Kaplan and Meier and another due to Nelson and Åalen. The following summarizes these estimators without adjustments needed for truncation. To compute these estimators, you can do the following :

• Use the problem to identify the data and the non-censored losses/events. Specifically, begin with a data set ( {x_1, ldots, x_n } ). Let (t_1 < ldots < t_k ) be the subset of ( {x_1, ldots, x_n } ) that correspond to unique non-censored times. For each event time (non-censored losses) (t_j), let (s_j) be the number of events.
• For each event time (non-censored losses) (t_j), determine the risk set, defined to be
  $$R_j=sum_{i=1}^n I(x_i geq t_j) $$
• The Kaplan-Meier estimator of the survival function (S(x) = 1 – F(x) ) is defined to be
  $$widehat{S}_{KM}(x)=prod_{j:t_j leq x}left(1-frac{s_j}{R_j}right)$$
• The Nelson-Åalen estimator is defined to be
  $$widehat{S}_{NA}(x)=e^{-hat{H}(x)}$$based on the following estimate of the cumulative hazard function
  $$hat{H}(x)=sum_{j:t_j leq x} frac{s_j}{R_j} .$$
• To determine the reliability of the Kaplan-Meier estimator, you may use Greenwood’s approximation to the variance, given as:
  $$widehat{Var}(widehat{S}_{KM}(x))=(widehat{S}_{KM}(x))^2 sum_{j:t_j leq x} frac{s_j}{R_j(R_j-s_j)} $$

---

For right-censored data (without truncation), the following gives an expression for a nonparametric estimator of the survival function due to Nelson and Åalen. To compute this estimator:

• Use the problem to identify the data and the non-censored losses/events. Specifically, begin with a data set ( {x_1, ldots, x_n } ). Let (t_1 < ldots < t_k ) be the subset of ( {x_1, ldots, x_n } ) that correspond to unique non-censored times. For each event time (non-censored losses) (t_j), let (s_j) be the number of events.
• For each event time (non-censored losses) (t_j), determine the risk set, defined to be (R_j=sum_{i=1}^n I(x_i geq t_j) ).
• The Nelson-Åalen estimator of the survival function (S(x) = 1 – F(x) ) is defined to be
  $$widehat{S}_{NA}(x)=e^{-hat{H}(x)}$$which based on the following estimate of the cumulative hazard function
  $$hat{H}(x)=sum_{j:t_j leq x} frac{s_j}{R_j} .$$
• Here is another set of notation that facilitates a connection to the parametric estimator. For each observation (x_i), define ( delta_i ) to be a binary variable that indicates whether the observation was censored. With this, we can write the risk set at (x_j) to be (R(x_j)=sum_{i=1}^n delta_i I(x_i geq x_j) ) and the cumulative hazard as $$hat{H}(x)=sum_{j:x_j leq x} frac{delta_j}{R(x_j)} .$$

For a parametric fit, suppose that we have a density (f) and distribution function (F) that depends on ( theta ). Define the log-likelihood to be
$$ L(theta) =sum_{i=1}^n left{ (1-delta_i) ln f(x_i) + delta_i ln left(1-F(x_i) right) right}.$$ The maximum likelihood estimator (hat{theta}) is defined to be that value of ( theta ) that maximizes ( L(theta) ). Often, this is done by taking a derivative of ( L(theta) ), setting it equal to zero, and solving for ( theta ).

For Bayesian estimation problems:

• Use the problem to identify the conditional model distribution of each observation given (theta), (f_{X|theta}(x|theta)). Then get the likelihood of the model distribution, which is typically the product of the model distributions for all observations,
  $$f_{X|theta}(x_1…x_n|theta)=prod_{i=1}^N f_{X|theta}(x_i|theta) .$$
• From the question identify the prior distribution of (theta), (pi(theta)).
• Now obtain the posterior distribution using the Bayes theorem;
  $$pi(theta|x_1…x_n)=frac{f_{X|theta}(x_1…x_n|theta)pi(theta)}{int f_{X|theta}(x_1…x_n|theta)pi(theta) dtheta}$$Note that the denominator of this expression ( f(x) = int f_{X|theta}(x_1…x_n|theta)pi(theta) dtheta) is the marginal distribution of (x) and does not involve ( theta ).
• With this, to calculate the posterior probability that (a lt theta lt b);
  $$pi(a lt theta lt b|x_1…x_n)=int_{a}^{b} pi(theta|x_1…x_n) dtheta$$
• For a new data (y) , the predictive distribution is;
  $$f(y|x_1…x_n)=int f(y|theta)pi(theta|x_1…x_n)dtheta .$$Here, (f(y|theta)) is obtained from the conditional model distribution

There are several ways that one can use the predictive distribution. For some problems, you may wish to calculate Bayesian prediction for a new data (y). Specifically, given the past observations (x_1…x_n), you can use
$$E(y|x_1…x_n)=int E(y|theta)pi(theta|x_1…x_n)dtheta. $$Here, (E(y|theta)) is the expected value of the conditional model distribution.

Method of moments and percentile matching are estimation techniques that provide alternatives to the method of maximum likelihood. Because they are non-recursive, they also can be used to provide starting values to begin a maximum likelihood recursion.

Under the method of moments, one matches the empirical (non-parametric) moments to moments specified by a parametric distribution. Assume that you have a distribution with (c) parameters (typically, (c=1,2,3,4)). Then,

• Determine the parameters in terms of the moments (mathrm{E}~ X^k = mu_k^{prime} ) for ( k=1, 2, ldots, c).
• Determine the sample (also known as the empirical or nonparametric) moments given by
  $$m_k =frac{1}{n}sum_{i=1}^{n} x^{k}, text{for~} k=1, 2, ldots,c$$
• Approximating (m_k approx mu_k^{prime}) for (k=1, 2, ldots,c), solve for the parameters. The resulting parameter estimates are method of moment estimators.

In the same way, under the percentile matching, one matches the empirical (non-parametric) percentiles or quantiles to those specified by a parametric distribution. As before, assume that you have a distribution with (c) parameters. Then,

• Determine the parameters in terms of the quantiles (pi_{q_k} = F^{-1}(q_k)) for selected proportions (q_k) for ( k=1, 2, ldots, c).
• Determine the corresponding sample quantiles. It is common to use the smoothed empirical percentile given by
  $$ hat{pi}_q = (1-h) X_{(j)} + h X_{(j+1)}$$
  where (j=[(n+1)q]) and, (h=(n+1)q-j), and (X_{(1)}, ldots, X_{(n)}) are the ordered values (the order statistics) corresponding to (X_1, ldots, X_n).
• Approximating (hat{pi}_{q_k} approx pi_{q_k}) for (k=1, 2, ldots,c), solve for the parameters. The resulting parameter estimates are percentile matching estimators.