Existence of optimal weighted least squares estimate for 3-parametric exponential model

François DUBEAU and Youness MIR

Département de mathématiques Université de Sherbrooke 2500 Boulevard de l’Université Sherbrooke (Qc), Canada, J1K 2R1

francois.dubeau@usherbrooke.ca youness.mir@usherbrooke.ca

Abstract : Given a set of data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m} we present sufficient conditions which guarantee the existence of the weighted least squares estimate for a 3parametric exponential function. To this end, we suggest a classification of the data based on their increasing or decreasing trend, and on their convexity or concavity form. We discuss a choice of the initial approximation and present a numerical example.

Key Words : weighted least squares, parametric exponential function, least squares estimate, existence problem, classification of data.

AMS Subject Classification : 65D10, 62J02.

1. Introduction

Very often in experimental research data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m} are obtained and we have to find a good model to fit them. In this paper we consider the 3-parametric exponential model

(1.1)  f(t; a, b, c) = a + be^{ct}

where a, b, and c are real parameters. These parameters will be estimated from the data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m} where the t_{i}’s represent the values of the independent variable, the f_{i}’s are the values of the measured function, and each ω_{i} is the weight associated to the data (t_{i}, f_{i}).

It is well known that such model are useful in different fields of applied research such that economy [14], biology [11], physics [15], chemistry [1], medicine [9], engineering [11, 13] and agriculture [4].

To estimate the parameters a, b, and c we use the weighted least squares method. We minimize the functional

(1.2)  F (a, b, c) = \frac{1}{2} ∑_{i=1}^{m} ω_{i}(f(t_{i}; a, b, c) − f_{i})^{2} = \frac{1}{2} ∑_{i=1}^{m} ω_{i}(a + be^{ct_{i}} − f_{i})^{2}

on ℝ^{3} and we find the weighted least squares estimate (WLSE) given by

(a^{∗}, b^{∗}, c^{∗}) = arg min_{(a,b,c)∈ℝ^{3}} F (a, b, c).

In the sequel we assume that m ≥ 3 and the t_{i}’s are strictly increasing

t_{1} < t_{2} < ⋯ < t_{m}.

The goal of this paper is to present sufficient conditions for the existence of the WLSE for the 3-parametric exponential model (1.1) to complement the papers [5] and [6]. It appends that the proof contained in these papers contains some mistakes and we present a way to get the result correctly. Our approach is based on orthogonal polynomials to obtain a rough classification of the data and then deduce sufficient conditions for the existence of the WLSE. Other approaches concerning the existence of the WLSE exist. In particular, using the concept of ”existence level”, some theorems about the existence of the WLSE are obtained in [2] and [3].

In Section 2 we briefly review the polynomial least squares fitting and suggest an elementary classification of the data based on the coefficients of the best polynomial with respect to an orthogonal basis of polynomials. In Section 3 we present sufficient conditions for the existence of the weighted least square estimate for the 3-parametric exponential model. In Section 4 we suggest two ways to initialize the numerical algorithm to find the optimal solution. A numerical example follows in Section 5.

2. Polynomial least squares fitting and data classification

Let P_{n} be the set of polynomials of degree ≤ n. In this section, we look for the best n-degree polynomial fitting to the data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m}. This polynomial is given by

p_{n}^{∗} = arg min_{p∈P_{n}} G(p)

where

G(p) = \frac{1}{2} ∑_{i=1}^{m} ω_{i}(p(t_{i}) − f_{i})^{2}.

This polynomial exists and is unique as long as n + 1 ≤ m. Indeed in this case

G(p) = \frac{1}{2}∥p − f∥^{2}

where ∥ · ∥ is the norm on P_{n} induced by the scalar product defined by

⟨p, q⟩ = ∑_{i=1}^{m} ω_{i}p(t_{i})q(t_{i})

for any pair of polynomials p and q in P_{n}. For the f_{i}’s we use the notation f_{i} = f(t_{i}) (i = 1, ..., m). It is well known that p_{n}^{∗} is characterized by the normal equations

⟨p_{n}^{∗} − f, p⟩ = 0

for all p ∈ P_{n}.

In this setting, we can find a sequence of orthogonal polynomials by applying the Gram-Schmidt orthogonalization process to the standard basis {1, t, t^{2}, ⋯, t^{n}} of P_{n}. These orthogonal polynomials are given by

q_{0}(t) = 1, q_{1}(t) = t − α_{1},

and for j ≥ 2

q_{j}(t) = (t − α_{j})q_{j−1}(t) − β_{j}q_{j−2}(t)

where

α_{j} = \frac{⟨tq_{j−1}, q_{j−1}⟩}{⟨q_{j−1}, q_{j−1}⟩} (j = 1, 2, 3, …),

and

β_{j} = \frac{⟨tq_{j−1}, q_{j−2}⟩}{⟨q_{j−2}, q_{j−2}⟩} (j = 2, 3, 4, …).

Hence the best n-degree least squares polynomial p_{n}^{∗} can be written as

p_{n}^{∗}(t) = ∑_{j=0}^{n} γ_{j}^{∗}q_{j}(t)

where

γ_{j}^{∗} = \frac{⟨f, q_{j}⟩}{⟨q_{j}, q_{j}⟩} (j = 0, 1, …, n).

Let us observe that

p_{n}^{∗}(t) = γ_{n}^{∗}q_{n}(t) + p_{n−1}^{∗}(t).

In the next section we will use the following result.

Lemma 2.1. ⟨p_{n−1}^{∗} − f, t^{n}⟩ = −γ_{n}^{∗}∥q_{n}∥^{2} for n = 1, …, m − 1 .

Proof. Since q_{n}(t) = t^{n} + q̃_{n−1}(t) where q̃_{n−1}(t) is a polynomial of degree ≤ n − 1, we have

γ_{n}^{∗}∥q_{n}∥^{2} = ⟨γ_{n}^{∗}q_{n}, q_{n}⟩
= ⟨p_{n}^{∗} − p_{n−1}^{∗}, q_{n}⟩
= ⟨p_{n}^{∗} − f, q_{n}⟩ + ⟨f − p_{n−1}^{∗}, q_{n}⟩
= −⟨p_{n−1}^{∗} − f, t^{n} + q̃_{n−1}⟩
= −⟨p_{n−1}^{∗} − f, t^{n}⟩.  □

Lemma 2.2. If the q_{j}’s are the orthogonal polynomials associated to {(ω_{i}, t_{i})}_{i=1}^{m}, the orthogonal polynomials q̃_{j}’s associated to {(ω_{i}, t̃_{i} = −t_{i})}_{i=1}^{m} are given by

q̃_{j}(t) = (−1)^{j}q_{j}(−t). □

Definition 2.3. The data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m} are said to be

(i) stationary if γ_{1}^{∗} = 0; (ii) essentially increasing, respectively decreasing, if γ_{1}^{∗} > 0, respectively γ_{1}^{∗} < 0; (iii) essentially linear, if γ_{2}^{∗} = 0; (iv) essentially convex, respectively concave, if γ_{2}^{∗} > 0, respectively γ_{2}^{∗} < 0. □

Remark 2.4. If the data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m} are essentially increasing, resp. decreasing, then the data {(ω_{i}, −t_{i}, f_{i})}_{i=1}^{m} are essentially decreasing, resp. increasing. Also if the data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m} are essentially convex, resp. concave, then the data {(ω_{i}, t_{i}, −f_{i})}_{i=1}^{m}are essentially concave, resp. convex. □

Let us define the following expression

f_{•j} = \frac{∑_{i=1,i≠j}^{m} ω_{i}f_{i}}{∑_{i=1,i≠j}^{m} ω_{i}}

for j = 1, ..., m.

Definition 2.5. The data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m} are said to be strongly regular if they are (i) essentially increasing and convex, and f_{m} > f_{•m} and f_{m−1} > f_{•m}, or (ii) essentially decreasing and convex, and f_{1} > f_{•1} and f_{2} > f_{•1}, or (iii) essentially increasing and concave, and f_{1} < f_{•1} and f_{2} < f_{•1}, or (iv) essentially decreasing and concave, and f_{m} < f_{•m} and f_{m−1} < f_{•m}. □

3. Existence of the least squares estimate

In this section we discuss the weighted least squares estimate existence problem. We give conditions under which there exists a (a^{∗}, b^{∗}, c^{∗}) ∈ ℝ^{3} such that

(a^{∗}, b^{∗}, c^{∗}) = arg min_{(a,b,c)∈ℝ^{3}} F (a, b, c).

It is necessary to add some conditions on the data for this existence problem as shown by the following example.

Example 3.1. Suppose the data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m} are such that the data points {(t_{i}, f_{i})}_{i=1}^{m} are all on a given line y = kt + l.

If k = 0 hence l = f_{i} for i = 1, ..., m, and let (a^{∗}, b^{∗}, c^{∗}) ∈ ℝ^{3} such that a^{∗} + b^{∗}e^{c^{∗}t_{i}} = l for i = 1, ..., m. Then

a^{∗} = l, b^{∗} = 0, and c^{∗} ∈ ℝ^{3},

or

a^{∗} + b^{∗} = l, and c^{∗} = 0,

are all optimal solutions such that F (a^{∗}, b^{∗}, c^{∗}) = 0.

If k ≠ 0, since the graphs of y = f(t; a, b, c) and y = kt + l intersect in at most two points and m ≥ 3, we have

F (a, b, c) = \frac{1}{2} ∑_{i=1}^{m} ω_{i}(a + be^{ct_{i}} − f_{i})^{2} > 0.

But for

(a, b, c) = (f_{0} − \frac{1}{c}, \frac{e^{−ckt_{0}}}{c} , ck) ∈ ℝ^{3}

where f_{0} = kt_{0} + l and c ≠ 0, we have

lim_{c→0} F (f_{0} − \frac{1}{c}, \frac{e^{−ckt_{0}}}{c} , ck) = lim_{c→0} \frac{1}{2} ∑_{i=1}^{m} ω_{i}(f_{0} + \frac{e^{ck(t_{i}−t_{0})} − 1}{c} − f_{i})^{2}

= \frac{1}{2} ∑_{i=1}^{m} ω_{i}(f_{0} + k(t_{i} − t_{0}) − f_{i})^{2}

= \frac{1}{2} ∑_{i=1}^{m} ω_{i}(kt_{i} + l − f_{i})^{2}
= 0.

Hence

inf_{(a,b,c)∈ℝ^{3}} F (a, b, c) = 0

and the problem has no solution. □

To obtain an existence result we will assume that the data satisfy certain sufficient conditions. For example the data will be essentially increasing or decreasing and essentially convex or concave.

Let us consider the following subsets of ℝ^{3}

U_{++} = { (a, b, c) ∈ ℝ^{3} | b ≥ 0 and c ≥ 0 },
U_{+−} = { (a, b, c) ∈ ℝ^{3} | b ≥ 0 and c ≤ 0 },
U_{−+} = { (a, b, c) ∈ ℝ^{3} | b ≤ 0 and c ≥ 0 },
U_{−−} = { (a, b, c) ∈ ℝ^{3} | b ≤ 0 and c ≤ 0 }.

The boundary of these subsets are

∂U_{++} = Γ_{b}^{+} ∪ Γ_{c}^{+},  ∂U_{+−} = Γ_{b}^{+} ∪ Γ_{c}^{−},
∂U_{−+} = Γ_{b}^{−} ∪ Γ_{c}^{+},  ∂U_{−−} = Γ_{b}^{−} ∪ Γ_{c}^{−},

where

Γ_{b}^{+} = { (a, b, 0) ∈ ℝ^{3} | b ≥ 0 },
Γ_{b}^{−} = { (a, b, 0) ∈ ℝ^{3} | b ≤ 0 },
Γ_{c}^{+} = { (a, 0, c) ∈ ℝ^{3} | c ≥ 0 },
Γ_{c}^{−} = { (a, 0, c) ∈ ℝ^{3} | c ≤ 0 }.

For the interior of a subset U of ℝ^{3}, defined by U \ ∂U, we will use the notation Int(U).

Theorem 3.2. Let the data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m} be essentially convex (γ_{2}^{∗} > 0) and strongly regular. (A) If the data are essentially increasing (γ_{1}^{∗} > 0) then there exists (a^{∗}, b^{∗}, c^{∗}) ∈ Int(U_{++}) such that

F (a^{∗}, b^{∗}, c^{∗}) = min_{(a,b,c)∈U_{++}} F (a, b, c).

(B) If the data are essentially decreasing (γ_{1}^{∗} < 0) then there exists (a^{∗}, b^{∗}, c^{∗}) ∈ Int(U_{+−}) such that

F (a^{∗}, b^{∗}, c^{∗}) = min_{(a,b,c)∈U_{+−}} F (a, b, c).

Proof. We assume that the t_{i}’s are as follows

t_{1} < t_{2} = 0 < t_{3} < ⋯ < t_{m}.

(A) Since F (a, b, c) ≥ 0 then

0 ≤ inf_{(a,b,c)∈U_{++}} F (a, b, c) = F^{∗} < +∞.

Step 1. Let us show that F^{∗} cannot be a value of F on the boundary ∂U_{++} = Γ_{b}^{+} ∪ Γ_{c}^{+}. On Γ_{b}^{+}, c = 0 and

F (a, b, 0) = \frac{1}{2} ∑_{i=1}^{m} ω_{i}(a + b − f_{i})^{2} ≥ \frac{1}{2} ∑_{i=1}^{m} ω_{i}(γ_{0}^{∗} − f_{i})^{2} = F_{0}.

This value can be obtained for any (ã, b̃, 0) ∈ Γ_{b}^{+} such that ã + b̃ = γ_{0}^{∗} and b̃ ≥ 0. Also, because γ_{1}^{∗} > 0, and

\frac{∂F}{∂c} (ã, b̃, 0)_{{|ã+b̃=γ_{0} \atop b̃>0}} = b̃ ∑_{i=1}^{m} ω_{i}(γ_{0}^{∗} − f_{i})t_{i}
= b̃⟨p_{0}^{∗} − f, t⟩
= −b̃γ_{1}^{∗}∥q_{1}∥^{2}
< 0

for b̃ > 0, there exists a point (ã, b̃, c̃) ∈ Int(U_{++}) such that

F (ã, b̃, c̃)|_{{ã+b̃=γ^{∗}_{0} \atop b̃>0}} < F (ã, b̃, 0)|_{{ã+b̃=γ^{∗}_{0} \atop b̃>0}} < F (a, b, 0)_{|(a,b,0)∈Γ_{b}^{+}}.

On Γ_{c}^{+}, b = 0 and

F (a, 0, c) = \frac{1}{2} ∑_{i=1}^{m} ω_{i}(a − f_{i})^{2}

≥ \frac{1}{2} ∑_{i=1}^{m} ω_{i}(γ_{0}^{∗} − f_{i})^{2}
≥ F (a, b, 0)|_{{ã+b̃=γ^{∗}_{0} \atop b̃>0}}
> F (ã, b̃, c̃)|_{{ã+b̃=γ^{∗}_{0} \atop b̃>0}} .

Hence F^{∗} cannot be attained on ∂U_{++}.

Step 2. Now, let {a_{n}, b_{n}, c_{n}}_{n=1}^{+∞} be a sequence in U_{++} such that

0 ≤ lim_{n→+∞} F (a_{n}, b_{n}, c_{n}) = inf_{(a,b,c)∈U_{++}} F (a, b, c) = F^{∗} < +∞.

Let

r_{i}(n) = a_{n} + b_{n}e^{c_{n}t_{i}} − f_{i}

for i = 1, ..., m.

Without loss of generality, whenever we have an unbounded sequence we may assume that it tends to +∞ or to −∞, because we can take a subsequence. Similarly, whenever we have a bounded sequence, we may assume it is convergent, otherwise by the Bolzano-Weierstrass theorem, we can take a convergent subsequence.

Let us show that the sequence {a_{n}, b_{n}, c_{n}}_{n=1}^{+∞} is bounded by showing that F (a_{n}, b_{n}, c_{n}) cannot attain its infimum F^{∗} in either the following two cases :

2.I |a_{n}| → +∞ ; 2.II lim_{n→+∞} a_{n} = a^{∗} and lim_{n→+∞}(b_{n} + c_{n}) = +∞ .

Case 2.I. Since

r_{2}(n) = a_{n} + b_{n} − f_{2}

we must have lim_{n→+∞} b_{n} = − lim_{n→+∞} a_{n} and, from the fact that b_{n} ≥ 0, it follows that lim_{n→+∞} a_{n} = −∞ and lim_{n→+∞} b_{n} = +∞. We can suppose that lim_{n→+∞}(a_{n} +b_{n}) = l ∈ ℝ if not lim_{n→+∞} r_{2}^{2}(n) = +∞. If the sequence {c_{n}}_{n=1}^{+∞} is unbounded, then lim_{n→+∞} c_{n} = +∞. Take any t_{i_{0}} in {t_{3}, t_{4}, ..., t_{m}}, then c_{n}t_{i_{0}} > 0. Since b_{n} > 0 we have

a_{n} + b_{n}e^{c_{n}t_{i_{0}}} ≥ a_{n} + b_{n}(1 + c_{n}t_{i_{0}}) = (a_{n} + b_{n}) + b_{n}c_{n}t_{i_{0}}

because e^{c_{n}t_{i_{0}}} ≥ 1 + c_{n}t_{i_{0}}, and

lim_{n→+∞}(a_{n} + b_{n}) + b_{n}c_{n}t_{i_{0}} = l + ∞ = +∞.

Hence lim_{n→+∞} c_{n} = c̃ ≥ 0. If c̃ > 0, as before we obtain lim_{n→+∞} a_{n} + b_{n}e^{c_{n}t_{i}} = +∞. If c̃ = 0, from the mean value theorem we have

a_{n} + b_{n}e^{c_{n}t_{i}} = (a_{n} + b_{n}) + b_{n}c_{n}t_{i}e^{θ_{i}(n)c_{n}t_{i}}

where θ_{i}(n) ∈ (0, 1) for all i = 1, 2, ..., m and n = 1, 2, .... If the sequence {b_{n}c_{n}}_{n=1}^{+∞} is unbounded, as before we obtain lim_{n→+∞} a_{n} + b_{n}e^{c_{n}t_{i}} = +∞. Hence {b_{n}c_{n}}_{n=1}^{+∞} is bounded and we can suppose that it converges to k ≥ 0. It follows that

lim_{n→+∞} a_{n} + b_{n}e^{c_{n}t_{i}} = l + kt_{i}

for i = 1, ..., m. Hence

lim_{n→+∞} F (a_{n}, b_{n}, c_{n}) = \frac{1}{2} ∑_{i=1}^{m} ω_{i}(l + kt_{i} − f_{i})^{2} = F_{0} ≥ \frac{1}{2} ∑_{i=1}^{m} ω_{i}(p_{1}^{∗}(t_{i}) − f_{i})^{2}

where p_{1}^{∗} is the best least squares polynomial of degree ≤ 1 for the data. We have

p_{1}^{∗}(t) = γ_{0}^{∗}q_{0}(t) + γ_{1}^{∗}q_{1}(t)
= γ_{0}^{∗} + γ_{1}^{∗}(t − α_{1})
= (γ_{0}^{∗} − γ_{1}^{∗}α_{1}) + γ_{1}^{∗}t
= l^{∗} + k^{∗}t

where l^{∗} = γ_{0}^{∗} − γ_{1}^{∗}α_{1} and k^{∗} = γ_{1}^{∗} > 0. Let us choose any value (τ_{0} and define a continuous function ϕ : ℝ ↦ ℝ by the formula

ϕ(λ) = \{ \begin{matrix}\frac{1}{2} ∑_{i=1}^{m} ω_{i}(k^{∗}τ_{0} + \frac{e^{λk^{∗}(t_{i}−τ_{0})} − 1}{λ} − (f_{i} − l^{∗}))^{2} & if & λ ≠ 0, \\ \frac{1}{2} ∑_{i=1}^{m} ω_{i}(l^{∗} + k^{∗}t_{i} − f_{i})^{2} & if & λ = 0. \end{matrix}

Let us observe that

ϕ(λ) = F (k^{∗}τ_{0} + l^{∗} − \frac{1}{λ}, \frac{e^{−λk^{∗}τ_{0}}}{λ} , λk^{∗})

and

(k^{∗}τ_{0} + l^{∗} − \frac{1}{λ}, \frac{e^{−λk^{∗}τ_{0}}}{λ} , λk^{∗}) ∈ Int(U_{++})

for λ > 0. Its derivative is given by

ϕ^{(1)}(λ) = \{ \begin{matrix}∑_{i=1}^{m} ω_{i}(k^{∗}τ_{0} + \frac{e^{λk^{∗}(t_{i}−τ_{0})} − 1}{λ} − (f_{i} − l^{∗}))(\frac{e^{λk^{∗}(t_{i}−τ_{0})}(λk^{∗}(t_{i} − τ_{0}) − 1) + 1}{λ^{2}}) \\ if λ ≠ 0, \\ \frac{(k^{∗})^{2}}{2} ∑_{i=1}^{m} ω_{i}(l^{∗} + k^{∗}t_{i} − f_{i})(t_{i} − τ_{0})^{2} if λ = 0 \end{matrix}

is also continuous on ℝ. But

ϕ^{(1)}(0) = \frac{(k^{∗})^{2}}{2} ⟨l^{∗} + k^{∗}t − f, (t − τ_{0})^{2}⟩ = \frac{(k^{∗})^{2}}{2} ⟨l^{∗} + k^{∗}t − f, t^{2}⟩ = −\frac{(k^{∗})^{2}}{2} γ_{2}^{∗}∥q_{2}∥^{2}.

Since γ_{2}^{∗} > 0, it follows that ϕ^{(1)}(0) < 0 and the function ϕ is decreasing on a neighborhood of λ = 0. Hence their exists a λ̃ > 0 such that ϕ(λ̃) < ϕ(0) ≤ F_{0} and it follows that F^{∗} < F_{0} and the sequence {a_{n}}_{n=1}^{+∞} is bounded.

Case 2.II. Now we show that the sequence {(b_{n}, c_{n})}_{n=1}^{+∞} is bounded. To show that we consider the 4 possible situations

(1) b_{n} → +∞ and c_{n} → +∞ ,
(2) b_{n} → +∞ and c_{n} → c^{∗} ≥ 0 ,
(3) b_{n} → b^{∗} > 0 and c_{n} → +∞ ,
(4) b_{n} → b^{∗} = 0 and c_{n} → +∞ .

For the cases (1), (2), and (3), we observe that {F (a_{n}, b_{n}, c_{n})}_{n=1}^{+∞} diverge to +∞ because lim_{n→+∞} r_{2}(n) = +∞. It remains to consider only the last case.

Let us assume that lim_{n→+∞} a_{n} = a^{∗}, lim_{n→+∞} b_{n} = b^{∗} = 0, and lim_{n→+∞} c_{n} = +∞. We also suppose that

lim_{n→+∞} b_{n}e^{c_{n}t_{m}} = L < +∞

because in the contrary lim_{n→+∞} r_{m}^{2}(n) = +∞. We have

lim_{n→+∞} b_{n}e^{c_{n}t_{i}} = lim_{n→+∞} b_{n}e^{c_{n}t_{m}} lim_{n→+∞} e^{c_{n}(t_{i}−t_{m})} = L · 0 = 0

for i = 1, ..., m − 1. Hence

F_{1} = lim_{n→+∞} F (a_{n}, b_{n}, c_{n})
= \frac{1}{2} ∑_{i=1}^{m−1} ω_{i}(a^{∗} − f_{i})^{2} + \frac{1}{2}ω_{m}(L − f_{m})^{2}

≥ \frac{1}{2} ∑_{i=1}^{m−1} ω_{i}(f_{•m} − f_{i})^{2}.

We look for a point (f_{•m}, b′, c′) in Int(U_{++}) such that

F (f_{•m}, b′, c′) < \frac{1}{2} ∑_{i=1}^{m−1} ω_{i}(f_{•m} − f_{i})^{2} ≤ F_{1}.

Using the strong monotony for f_{m}, we consider the curve sketched in Figure 1 and defined by

Γ = { (f_{•m}, (f_{m} − f_{•m})e^{−ct_{m}}, c) ∈ U_{++} | c ≥ 0}.

Figure 1. The curve Γ.

Let

ψ(c) = F (f_{•m}, (f_{m} − f_{•m})e^{−ct_{m}}, c) = \frac{1}{2} ∑_{i=1}^{m} ω_{i}(f_{•m} + (f_{m} − f_{•m})e^{c(t_{i}−t_{m})} − f_{i})^{2}.

Then

lim_{c→+∞} ψ(c) = \frac{1}{2} ∑_{i=1}^{m−1} ω_{i}(f_{•m} − f_{i})^{2}.

Also

ψ^{(1)}(c) = ∑_{i=1}^{m−1} ω_{i}(f_{•m} + (f_{m} − f_{•m})e^{c(t_{i}−t_{m})} − f_{i})(f_{m} − f_{•m})(t_{i} − t_{m})e^{c(t_{i}−t_{m})},

and if we consider the ratio

\frac{ψ^{(1)}(c)}{e^{c(t_{m−1}−t_{m})}} = ∑_{i=1}^{m−1} ω_{i}(f_{•m}+(f_{m}−f_{•m})e^{c(t_{i}−t_{m})}−f_{i})(f_{m}−f_{•m})(t_{i}−t_{m−1})e^{c(t_{i}−t_{m−1})},

we have

lim_{c→+∞} \frac{ψ^{(1)}(c)}{e^{c(t_{m−1}−t_{m})}} = ω_{m−1}(t_{m} − t_{m−1})(f_{m} − f_{•m})(f_{m−1} − f_{•m}) > 0

because the data are strongly monotonic. Hence ψ^{(1)}(c) > 0 for large enough c and there exists a c̄ such that ψ(c) is strictly increasing on (c̄, +∞). Then there exists a c̃ ∈ (c̄, +∞) such that

ψ(c̃) = F (f_{•m}, (f_{m} − f_{•m})e^{−c̃t_{m}}, c̃) < lim_{c→+∞} ψ(c) ≤ F_{1}.

(B) From (A) we have a solution for the data {(ω_{i}, −t_{i}, f_{i})}_{i=1}^{m}. Then we only have to replace t by −t in the model, which changes c by −c. □

We can also prove the following theorem for essentially concave data.

Theorem 3.3. Let the data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{m} be essentially concave (γ_{2}^{∗} < 0) and strongly regular. (A) If the data are essentially increasing (γ_{1}^{∗} > 0) then there exists (a^{∗}, b^{∗}, c^{∗}) ∈ Int(U_{−−}) such that

F (a^{∗}, b^{∗}, c^{∗}) = min_{(a,b,c)∈U_{−−}} F (a, b, c).

(B) If the data are essentially decreasing (γ_{1}^{∗} < 0) then there exists (a^{∗}, b^{∗}, c^{∗}) ∈ Int(U_{−+}) such that

F (a^{∗}, b^{∗}, c^{∗}) = min_{(a,b,c)∈U_{−+}} F (a, b, c). □

4. Choice of the initial approximation

Numerical methods for minimizing the functional (1.2), which is a nonlinear weighted least squares problem, require an initial approximation (a_{0}, b_{0}, c_{0}) ∈ U which is as good as possible. A bad choice of the initial approximation can have bad consequences as a high number of iterations, a convergence to a local minimum or no convergence at all. A good initial approximation can leads to convergence in few iterations to the global minimum.

Here we compare two methods to find the initial approximation, both based on the choice of three points that represent the data and the model. We select three values (τ_{k}, φ_{k}) (for k = 1, 2, 3) and we solve the system of three nonlinear equations

(4.1)  φ_{k} = a + be^{cτ_{k}}  (k = 1, 2, 3).

Under the condition that τ_{2} = \frac{τ_{1}+τ_{3}}{2} we get

a = φ_{1} − \frac{(φ_{2} − φ_{1})^{2}}{φ_{3} − 2φ_{2} + φ_{1}} ,
b = \frac{(φ_{2} − φ_{1})^{2}}{φ_{3} − 2φ_{2} + φ_{1}} (\frac{φ_{2} − φ_{1}}{φ_{3} − φ_{2}} )^{2τ_{1}/(τ_{3}−τ_{1})},
c = \frac{2}{τ_{3} − τ_{1}} ln (\frac{φ_{3} − φ_{2}}{φ_{2} − φ_{1}} ).

1st method. Let ω = Σ_{i=1}^{m}ω_{i} and set 0 < α < 1/2. We first choose j_{1} and i_{3} such that

1 = i_{1} < j_{1}  and  ∑_{i=i_{1}}^{j_{1}} ω_{i} ≈ αω,

i_{3} < j_{3} = 1  and  ∑_{i=i_{3}}^{j_{3}} ω_{i} ≈ αω,

and we set

τ_{1} = \overline{t}_{1} = \frac{∑_{i=i_{1}}^{j_{1}} ω_{i}t_{i}}{∑_{i=i_{1}}^{j_{1}} ω_{i}}  and φ_{1} = \overline{f}_{1} = \frac{∑_{i=i_{1}}^{j_{1}} ω_{i}f_{i}}{∑_{i=i_{1}}^{j_{1}} ω_{i}}

τ_{3} = \overline{t}_{3} = \frac{∑_{i=i_{3}}^{j_{3}} ω_{i}t_{i}}{∑_{i=i_{3}}^{j_{3}} ω_{i}}  and φ_{3} = \overline{f}_{3} = \frac{∑_{i=i_{3}}^{j_{3}} ω_{i}f_{i}}{∑_{i=i_{3}}^{j_{3}} ω_{i}} .

Then we set τ_{2} = \frac{τ_{1}+τ_{3}}{2} and we determine 1 ≤ i_{2} ≤ j_{2} ≤ m such that

τ_{2} ≈ \overline{t}_{2} = \frac{∑_{i=i_{2}}^{j_{2}} ω_{i}t_{i}}{∑_{i=i_{2}}^{j_{2}} ω_{i}}  and  ∑_{i=i_{2}}^{j_{2}} ω_{i} ≈ αω

and we set

φ_{2} = \overline{f}_{2} = \frac{∑_{i=i_{2}}^{j_{2}} ω_{i}f_{i}}{∑_{i=i_{2}}^{j_{2}} ω_{i}}

We note by (ā_{0}, b̄_{0}, c̄_{0}) the value obtained in solving (4.1).

2nd method. Let

p_{2}^{∗}(t) = γ_{0}^{∗}q_{0}(t) + γ_{1}^{∗}q_{1}(t) + γ_{2}^{∗}q_{2}(t)

be the best least squares polynomial approximation to the data. The idea is to use only the part of the parabola that best represents the data. This part of the parabola is on the left side or on the right side of the minimum or maximum point of the parabola.

If γ_{1}^{∗}γ_{2}^{∗} < 0, the data are essentially increasing and concave or essentially decreasing and convex (hence the left side is appropriate here), we set

t_{max} = min\{t_{m}, \frac{α_{1} + α_{2}}{2} − \frac{γ_{1}^{∗}}{2γ_{2}^{∗}} \},

and

L_{T} = t_{m} − t_{1}
L_{0} = t_{m} − t_{max},
L_{1} = t_{max} − t_{1},

then

τ_{1} = t_{1} + λ(0.5 + ρ\frac{L_{0}}{L_{T}} )L_{1},
τ_{3} = t_{max} − λ(0.5 − ρ\frac{L_{0}}{L_{T}} )L_{1},
τ_{2} = \frac{τ_{1} + τ_{3}}{2} = \frac{t_{1} + t_{max}}{2} + λρ\frac{L_{0}L_{1}}{L_{T}}

If γ_{1}^{∗}γ_{2}^{∗} > 0, the data are essentially increasing and convex or essentially decreasing and concave (hence the right side is appropriate here), we set

t_{min} = max\{t_{1}, \frac{α_{1} + α_{2}}{2} − \frac{γ_{1}^{∗}}{2γ_{2}^{∗}} \},

and

L_{T} = t_{m} − t_{1}
L_{0} = t_{min} − t_{1},
L_{1} = t_{m} − t_{min},

then

τ_{1} = t_{min} + λ(0.5 − ρ\frac{L_{0}}{L_{T}} )L_{1},
τ_{3} = t_{m} − λ(0.5 + ρ\frac{L_{0}}{L_{T}} )L_{1},
τ_{2} = \frac{τ_{1} + τ_{3}}{2} = \frac{t_{1} + t_{max}}{2} − λρ\frac{L_{0}L_{1}}{L_{T}} .

Finally, we set φ_{k} = p_{2}^{∗}(τ_{k}) for k = 1, 2, 3. The parameters λ and ρ are such that 0 ≤ λ < 1, ρ ≥ 0 and (0.5−ρ\frac{L_{0}}{L_{T}} )L_{1} ≥ 0. We note by (â_{0}, b̂_{0}, ĉ_{0}) the value obtained in solving (4.1).

5. A numerical example

In this example we use the 3-parametric exponential model to fit m = 100 data points generated for i = 1, ..., m by

t_{i} = −10 + \frac{i}{5},
f_{i} = −5 + 4e^{0.15t_{i}} + r_{i},
r_{i} ∼ N (0, σ)  for σ = 2,

and the weights ω_{i} are given by the following assignment

ω_{i} = \{ \begin{matrix}1 & if & 0 ≤ u_{i} < 0.30, \\ 2 & if & 0.30 ≤ u_{i} < 0.55, \\ 3 & if & 0.55 ≤ u_{i} < 0.75, \\ 4 & if & 0.75 ≤ u_{i} < 0.90, \\ 5 & if & 0.90 ≤ u_{i} ≤ 1.00, \end{matrix}

where u_{i} ∼ U[0, 1] is a uniformly distributed random variable over [0, 1] and r_{i} ∼ N (0, σ) is a normal random variable with mean μ = 0 and standard deviation σ. The Levenberg-Marquardt’s method [7], [8] has been used to find the optimal solution.

Model Equation optimal value

Optimal solution

Strait line  y = γ_{1}^{∗}q_{1} + γ_{0}^{∗}q_{0}  1109.08  \begin{matrix}γ_{0}^{∗} = 0.19892 \\ γ_{1}^{∗} = 0.62588 \end{matrix}

Parabola  y = γ_{2}^{∗}q_{2} + γ_{1}^{∗}q_{1} + γ_{0}^{∗}q_{0} 816.32  \begin{matrix}γ_{0}^{∗} = 0.19892 \\ γ_{1}^{∗} = 0.62588 \\ γ_{2}^{∗} = 0.05245 \end{matrix}

Exponential  y = a^{∗} + b^{∗}e^{c^{∗}t}  758.95  \begin{matrix}a^{∗} = −3.58159 \\ b^{∗} = 2.17507 \\ c^{∗} = 0.19784 \end{matrix}

Table 1. Optimal values of the best polynomial approximation of degree 1 and 2 and the exponential function y = a + be^{ct}.

To check the hypothesis of Theorem 3.2 and 3.3, Table 1 contains the coefficients of p_{2}^{∗}. For this example, the data are essentially increasing γ_{1}^{∗} = 0.62588 > 0 and

(a) Data  (b) Strait line

(c) Parabola  (d) Exponential function

Figure 2. Graphic representations of the models of Table 1 with the data {(ω_{i}, t_{i}, f_{i})}_{i=1}^{i=m}.

essentially convex γ_{2}^{∗} = 0.05245 > 0, hence the solution will be in Int(U_{++}). This table also contains the optimal value of the criteria \frac{1}{2}∥p_{i} − f∥^{2}for polynomials of degree i = 1 and i = 2. It also contains the optimal value and the optimal solution for the 3-parameters exponential model. Figure 2 contains the graphs of the data and the three models : strait line, parabola, and 3-parametric exponential model.

Initial approximation

\begin{matrix}ā_{0} = −3.72764 \\ b̄_{0} = 2.06578 \\ c̄_{0} = 0.21350 \end{matrix} \begin{matrix}â_{0} = −4.94429 \\ b̂_{0} = 3.77342 \\ ĉ_{0} = 0.13789 \end{matrix}

F (a_{0}, b_{0}, c_{0})  780.66  792.30
number of iterations  5  7
execution time (sec)  0.28  0.33
∇F (x^{∗})  ( \begin{matrix}1.283D − 07 \\ 5.271D − 07 \\ 1.570D − 05 \end{matrix} ) ( \begin{matrix}2.462D − 08 \\ −4.255D − 07 \\ −7.609D − 06 \end{matrix} )

Stopping criteria

∥∇F (x^{∗})∥_{∞} ≤ 0.5 · 10^{−4} 1.570D-05 7.609D-06

Table 2. Performance comparison of the method with respect to the two initial approximations.

In Table 2, we compare the performance of the method with respect to the choice of the initial solution determined by the two methods described in the preceding section. We have indicated the value of the criteria for the starting point, the number of iteration needed by the Levenberg-Marquardt’s method to obtain the optimal solution using this initial approximation, and the time required to obtain the optimal solution. The sup norm of the gradient at the optimum and the stopping criteria. Let us observe that for the first choice with α = 1/5, we get αω = 52.4, j_{1} = 21, i_{3} = 77, i_{2} = 41 et j_{2} = 60. For the second choice the results given are for λ = 3/8 and ρ = 2.

6. Acknowledgments

This work has been supported by a NSERC (Natural Sciences and Engineering Research Council of Canada) individual discovery grant of the first author.

References

[1] H.G. Bock, Randwetproblemmethoden zur parameteridentifizierung in systemen nichtlinearer differentialgleichungen, Ph.D. Thesis, Bonn, 1985.

[2] E.Z. Demidenko, On the existence of the least squares estimate in nonlinear growth curve models of exponential type, Comm. Statist. Theory Methods, 25 (1996), 159-182.

[3] E.Z. Demidenko, Is this the least squares estimates ?, Biometrika, 87 (2000), 437-452.

[4] R.C. Jain, R. Agrawal, and K.N. Singh, A whitin year growth model for crop yield forecasting, Biomed. J., 34 (1992), 501-511.

[5] D. Jukić and R. Scitovski, The best least squares approximation problem for a 3-parametric exponential regression model, ANZIAM J., 42 (2000), 254-266.

[6] D. Jukić, A necessary and sufficient criteria for the existence of the least squares estimate for a 3-parametric exponential function, Applied Mathematics and Computation, 147 (2004), 1-17.

[7] K. Levenberg, A method for the solution of certain nonlinear problems in least-squares, Quart. Appl. Math., 2 (1944), 164-166.

[8] D.W. Marquart, An algorithm for least-squares estimation of nonlinear inequalities, SIAM. J. Appl. Math., 11 (1963), 431-441.

[9] M. Marus̆ić and Z. Bajzer, Generalized two-parameter equation of growth, J. Math. Anal. Appl., 179 (1993), 446-462.

[10] H. Mühlig, Lösung praktischer approximationahfgaben durch parameteridentifikation, ZAMM, 73 (1993), T837-T839.

[11] J.D. Murray, Mathematical Biology, Springer, Berlin, 1989.

[12] R. Scitovski and D. Jukić, Total least squares problem for exponential function, Inverse Problems, 12 (1996), 341-349.

[13] R. Scitovski and D. Jukić, A method for solving the parameter identification problem for ordinary differential equations of the second order, Appl. Math. Comp., 74 (1996), 273-291.

[14] R. Scitovski and S. Kosanović, Rate of change in economics research, Economics Analysis, 19 (1985), 65-73.

[15] H. Späth, Numerik Vieweg, Braunschweig, 1994.