4.7. Appendices

The method of images is an intuitively appealing approach to boundary crossing problems. It involves the placement of a source term on one side of a boundary and a sink term on the other side (Daniels, 1982). The sink term has the effect of drawing off density from the source term as it reaches the boundary. The approach produces some nice generalities about boundary crossing problems.

With natural boundary conditions, it was shown previously in this chapter that the solution to the advection-diffusion equation with constant coefficients, , is

(4.23)

with initial conditions . An absorbing boundary at x = L can be achieved by placing a sink term with weight at x = 2L (Daniels, 1982). f(x, t) can then be expressed as

. (4.24)

Since f(x, t) vanishes at x = L, can be solved for by setting f(L, t) = 0 in equation (4.24). is obtained by completing the square in the second term.

appendix 4.b

In this appendix, I will discuss a method involving Laplace transforms for determining first passage distributions (Riccardia, 1977). I will develop the general approach and then show it can be used in the case of the Wiener drift process with a simple boundary to produce equation (4.7). I have based the derivation on several references: Darling and Siegert (1953), Prabhu (1965), Riccardia (1977), and Chhikara and Folks (1989).

This method takes advantage of the following theorem (due to Siegert (1951)). First, let X(t) be a homogeneous Markov process with continuous sample paths. Define f (x, t | x₀, t₀) as the conditional probability density function for X(t) = x given that X (t₀) = x₀. Also, define g (t, L | t₀, x₀) as the probability density function for the time T when X first reaches the state L, L > x₀. The random variable T can be expressed as:

. (4.25)

Also, let f^* denote the Laplace transform of f. In other words,

. (4.26)

Similarly, let g^* be the Laplace transform of g.

The following theorem Siegert (1951) is useful in determining g from f and L.

Theorem 4.1

If x₀ < L < x, then

. (4.27)

proof:

The proof follows by considering paths that lie at x > L at time t. Paths that are beyond L at time t must have first reached L at some time s with s < t. Thus, we can write the conditional probability distribution for x in terms of possible paths from x₀ to x:

. (4.28)

Applying the convolution theorem for Laplace transforms to equation (4.28) yields:

. (4.29)

Rearranging terms results in equation (4.27).

This theorem is useful in cases where the Laplace transform is known for f. The Laplace transform for g can then be determined from equation (4.27). If the inverse Laplace transform is known for g^*, then g can be obtained. In many cases, this is not practical. In some cases, such as the Wiener process and Weiner process with drift, the pertinent Laplace transforms and inverse Laplace transforms are known, and equation (4.27) can be used to determine the first passage distributions.

In the case of the Wiener Process with drift,

. (4.30)

After combining the exponents, completing the square, and rearranging terms, we end up with:

. (4.31)

Integrating (4.31) by parts yields:

. (4.32)

Substituting L for x₀ in equation (4.32) and plugging this into equation (4.27) yields:

. (4.33)

Making use of the fact that (Haberman, 1987), we arrive at:

. (4.34)

An alternative approach for determining the Laplace transform of the arrival time distribution is as follows. This approach begins by considering equation (4.5), the probability density for X(t) = x given that the process hasn't reached the barrier by time t. This can be written as for x < L, t < T (recall from above that T is defined as the time of absorption at the boundary). Also define

, (4.35)

and note that = prob(). P satisfies the backward Chapman-Kolmogorov equation (Cox and Miller, 1965); in other words, P satisfies

. (4.36)

Also, P can be related to the arrival time distribution by the relation

. (4.37)

Plugging equation (4.37) into equation (4.36) and taking the Laplace transform of both sides yields

. (4.38)

This is a second order linear ordinary differential equation and can be solved directly. The general form of the solution is

. (4.39)

₁and ₂ are the roots of the characteristic equation . Thus

,

and ₁ is positive (and real) and ₂ is negative. The particular solution can be obtained from the following information. First, g^* is bounded for > 0:

. (4.40)

Thus the coefficient B = 0, or else the second term of the general solution would become unbounded as . The coefficient A can be determined by noting that when x₀ = L, absorption is immediate, and g^* (t | ) = 1. This yields A = e^-L, and

. (4.41)

This is the same as equation (4.33) and is inverted in the same manner to give g(t).

appendix 4.c

In computing equation (4.8) an exponential overflow problem can be encountered. This equation involves multiplying an exponential that is large (sometimes larger than the machine can handle) by a standard normal probability that is very small. Dennis et al. (1991) present a method for combining these two terms in a numerical approximation for equation (4.8). Since I use this approximation extensively in computations, I will present the details.

The cumulative distribution function for the inverse Gaussian is

, (4.42)

where is the cdf of the standard normal distribution. Problems may arise in evaluating the above equation because the second term involves multiplying a large number, exp(·), by a very small number, (·). For certain combinations of r, , L and t either of these numbers may be beyond the precision of the computer used. Dennis, et al. (1991) present a method that circumvents this problem by combining the two components of the second term. The following has been modified from their approach.

First, making use of the relation , rewrite equation (4.42) as

, (4.43)

with

.

Now denote the pdf of the standard normal distribution as . It is easy to show that

, (4.44)

and thus

. (4.45)

It follows that

. (4.46)

This can be evaluated using the following approximations for (Abramowitz and Stegun, 1965). If x < 4,

, (4.47)

where q_y = 1/(1+d₀), d₀ = 0.2316419, d₁ = 0.319381530, d₂ = -0.356563782, d₃ = 1.781477937, d₄ = -1.821255978, and d₅ = 1.330274429. For

, . (4.48)

Therefore, for z < 4, using equations (4.46) and (4.47),

. (4.49)

And for , using equations (4.46) and (4.48),

. (4.50)

The code to evaluate G(t) using the above procedure is contained in Appendix 3.

appendix 4.d

Several of the simulations I perform require the generation of inverse Gaussian random variates. A standard procedure for generating random variates, x, from a probability density function, f(x), is to use the inverse of the cumulative distribution function, F(x). The procedure involves generating a uniform random variate on the range (0,1) and using the transformation x = F^-1(x) to generate the random variate. To perform this a closed form solution of F^-1 is required, but this is not known for the inverse Gaussian distribution.

Another approach utilizes a known relationship v = g(x), with the random variate, v, coming from an easily generated distribution. This procedure becomes a bit more complicated when there is more than one root, x_i, for a given observation, v₀, and it must be determined how to chose one of the roots.

Michael, et al. (1976) present a procedure for generating inverse Gaussian random variates using a transformation that yields two roots. One of the roots is selected with an assigned probability. The inverse Gaussian can be written as

, . (4.51)

The parameter , and . With the transformation

, (4.52)

V is distributed as ² with one degree of freedom. The ² variates, v, are easily generated as squares of standard normal random variates. For a particular observation, v₀, equation (4.52) has two roots:

. (4.53)

Michael, et al. (1976) show that the root x₁ should be chosen with probability

. (4.54)

Thus a general procedure for generating inverse Gaussian random variates is as follows. First, generate a random variate from a standard normal distribution and square it to generate an observation from ²₍₁₎, v₀. Next, use equation (4.53) to calculate the roots x₁ and x₂. Finally, perform a Bernoulli trial with equation (4.54) to select the appropriate root. The computer code to perform this procedure is provided in Appendix 3.

appendix 4.e

The results of the application of the two parameter travel time model (equation (4.7)) to continuous PIT tag data are contained in Table 4.4 - Table 4.6.