On Parameter Estimation by Nonlinear Least Squares in Some Special Two-Parameter Exponential Type Models

Download as pdf or txt
Download as pdf or txt
You are on page 1of 7

Appl. Math. Inf. Sci. 9, No.

6, 2925-2931 (2015) 2925


Applied Mathematics & Information Sciences
An International Journal

http://dx.doi.org/10.12785/amis/090620

On Parameter Estimation by Nonlinear Least Squares in


Some Special Two-Parameter Exponential Type Models
Darija Marković∗ and Luka Borozan
Department of Mathematics, J.J. Strossmayer University of Osijek, Trg Ljudevita Gaja 6, HR-31 000 Osijek, Croatia

Received: 20 Feb. 2015, Revised: 20 Apr. 2015, Accepted: 21 Apr. 2015


Published online: 1 Nov. 2015

Abstract: Two-parameter growth models of exponential type f (t; a, b) = g(t)exp(a + bh(t)), where a and b are unknown parameters
and g and h are some known functions, are frequently employed in many different areas such as biology, finance, statistic, medicine,
ect. The unknown parameters must be estimated from the data (wi ,ti , yi ), i = 1, . . . , n, where ti denote the values of the independent
variable, yi are respective estimates of regression function f and wi > 0 are some data weights. A very popular and widely used
method for parameter estimation is the method of least squares. In practice, to avoid using nonlinear regression, this kind of problems
are commonly transformed to linear, which is not statistically justified. In this paper we show that for strictly positive g and strictly
monotone h original nonlinear problem has a solution. Generalization in the l p norm (1 ≤ p < ∞) and some illustrative examples are
also given.

Keywords: two-parameter models, least squares, parameter estimation, existence problem, data fitting

1 Introduction Given models range from generally applied to very


specifically used.
In this paper we will investigate parameter estimation
problem for the models of the type:
2.1 Exponential regression
f (t; a, b) = g(t)ea+bh(t) , (1)
where a and b are unknown parameters and g and h are If we assume that the average rate of change of the
some known functions. This type of models are often population (P) over an interval of time is proportional to
used in applied research, such as biology, ecology, the size of the population (see [9,18]), we have the
political science, psychology, economics and finance (see following differential equation model:
e.g. [10,16,23,5,25]).
dP
The structure of the paper is as follows. In Section 2 = kP, (2)
we briefly describe few models of type (1). In Section 3 dt
we formulate ordinary least squares (OLS) fitting where (for growth) k is a positive constant.
problem for this type of models. In Section 4 we present To solve it, we can rewrite equation (2) and get a
our main result (Theorem 1) which guarantees the following equation
existence of the least squares estimate (LSE), provided
the data satisfy natural conditions. Illustrative numerical dP
examples are given in Section 5. = kdt.
P
Integration of both sides of the last equation yields
2 Some useful models of type (1) ln P = kt + C,

Now we will give brief descriptions of some models of or


type (1), which are commonly used in applied research. P(t) = ekt+C .
∗ Corresponding author e-mail: [email protected]
c 2015 NSP
Natural Sciences Publishing Cor.
2926 D. Marković, L. Borozan: On Parameter Estimation by Nonlinear Least Squares...

P (t) ✻ 2.3 Fox surplus-yield model

This model is mainly used in fishery sciences for


estimating the maximum sustainable yield. It was
proposed by Fox in 1970. (see [8]). The equilibrium
harvest or sustainable yield occurs when a fish stock’s
harvest rate H equals its growth rate G (see [6]), and that
is when
dN
✲ = G(N) − H(E, N) = 0. (3)
dt
t
The Fox model assumes a Gompertz growth of the
Fig. 1: Plots of the exponential model for some values of k and underlying stock size in the absence of harvest
C
K 
G(N) = rN ln ,
N
This model is sometimes also called log-level or where r is the intrinsic growth rate and K is the carrying
log-linear model (or regresion), since when dependent capacity of the environment. Another assumption of the
variable is log transformed relation becomes linear. Fox model is that harvest rate is proportional to fishing
Depending on the sign of parameter k, it can be increasing effort and the biomass of the stock; that is
or decreasing. Figure 1 gives some examples of the
exponential model.
H(E, N) = qEN,

where q is the catchability coefficient. Then, from (3), the


2.2 Power regression equilibrium biomass N ⋆ satisfies the relation
K
Assumption that percentage increases in independent
rN ⋆ ln − qEN ⋆ = 0
variable t leads to constant percentage changes in N⋆
dependant variable y (see e.g. [19,16]) implies:
which implies that the nonzero equilibrium is given by
dy dt
=k , q
y t N ⋆ (E) = Ke− r E .
and a simple calculation gives Therefore, the sustainable yield is
k lnt+C
y(t) = e . q
Y (E) = qEN ⋆ = qKEe− r E = EekE+C . (4)
This model has a property of scale invariance (see e.g.
[5,10,25]). Model is also known as power law, log-log The highest possible sustainable yield is called the
regression and allometric equation. maximum sustainable yield, denoted YMSY . The maximum
sustainable yield occurs when

dY q  q
y(t) ✻ = qK 1 − E e− r E = 0.
dE r

The corresponding optimal level of effort is given by


r
EMSY = .
q

Substituting this value of effort in (4), the maximum


sustainable yield is

t eC−1
YMSY = Kre−1 = − .
k
Fig. 2: Plots of the power model for some values of k and C

c 2015 NSP
Natural Sciences Publishing Cor.
Appl. Math. Inf. Sci. 9, No. 6, 2925-2931 (2015) / www.naturalspublishing.com/Journals.asp 2927

Y (E) ✻ data weights which describe the assumed relative


accuracy of the data.
There is no unique way to estimate the unknown
parameters in a nonlinear regression function and many
different methods have been proposed in literature (see
e.g. [1,2,3,11,12,20,22]).
If the errors in the measurements of the independent
variable are negligible, and the errors in the
measurements of the dependent variable are independent
✲ random variables following the normal distribution with
E
expectancy zero, i.e. that

Fig. 3: Plots of the Fox model for some values of k and C yi = f (ti ; a, b) + εi , i = 1, . . . , n,
then in practical applications the unknown parameters a
and b of model (1) are usually estimated in the sense of the
least squares (LS) method by minimizing the functional
2.4 Schumacher equation
n n
F2 (a, b) = ∑ wi εi2 = ∑ wi f (ti ; a, b) − yi
 2
Next model was proposed independently by Terazaki
i=1 i=1
[24], Johnson [14], Schumacher [21] and Michailoff [17].
The model assumes that the relative growth rate increases on the set P (R × (−∞, 0], R × [0, −∞) or R2 . A point
linearly with the squared inverse of time (see [4]): (a⋆ , b⋆ ) ∈ P such that F2 (a⋆ , b⋆ ) = inf(a,b)∈P F2 (a, b) is
called the least squares estimate (LS estimate), if it exists.
1 dy 1
=k 2. Special numerical methods have been developed for
y dt t the purpose of solving nonlinear LS problems (see e.g.
[7]). However, prior to minimization itself difficult
This leads us to model
questions are posed referring to the existence and
1
y(t) = eC−k t . uniqueness of the LS estimate as well as the problem of
determining a good initial approximation. In the next
This model is primarily employed for timber growth section we will prove existence result.
and yield modeling (see e.g. [4,13]). Models od type (1) belong to the class of models that
can be linearized by transforming the dependent
variables. This kind of transformation changes the error
structure, as well as the influences of the data values.
3 Formularization of the problem Because of that, it is not clear in which sense resulting
parameters would be optimal. However, result obtained
Parameters a and b of models of type (1) have to be from transformed data usually gives a good initial
estimated from the experimental or empirical data approximation for iterative minimization methods.
(wi ,ti , yi ), i = 1, . . . , n, n ≥ 3, where t1 < t2 < · · · < tn Models od type (1) also belong to the family of the
denote the values of the independent variable, yi are the quasilinear regression models. Recent result on parameter
respective measured function values and wi > 0 are the estimation problem for quasilinear models can be found in
[15].

y(t) ✻
4 The existence result
In this section we consider the existence of the best l p -
norm (1 ≤ p < ∞) estimator in regression model of the
form (1) where a and b are the unknown parameters, with
n
F(a, b) = ∑ wi g(ti )ea+bh(ti ) − yi
p
(5)
i=1

✲ The next lemma will be used in the proof of Theorem


t 1.
Lemma 1. Suppose we are given data (wi ,ti , yi ),
Fig. 4: Plots of the Schumacher model for some values of C and i = 1, . . . , n, n ≥ 3, such that t1 < t2 < . . . < tn and
k wi , yi > 0, i = 1, . . . , n.

c 2015 NSP
Natural Sciences Publishing Cor.
2928 D. Marković, L. Borozan: On Parameter Estimation by Nonlinear Least Squares...

A If h is strictly increasing, then: such that t1 < t2 < . . . < tn and wi , yi > 0, i = 1, . . . , n, then
(i) There exists a point in R × (−∞, 0) at which there exists a point (a⋆ , b⋆ ) ∈ P such that
functional F defined by (5) attains a value less
than ∑ni=2 wi yip . F(a⋆ , b⋆ ) = inf F(a, b).
(a,b)∈P
(ii) There exists a point in R × (0, ∞) at which
functional F defined by (5) attains a value less
p Proof. The proof will be carried out for the case when h is
than ∑n−1
i=1 wi yi .
B If h is strictly decreasing, then: strictly increasing. The proof for the second case (h is
(i) There exists a point in R × (0, ∞) at which strictly decreasing) is essentially identical, so it will be
functional F defined by (5) attains a value less omitted here.
than ∑ni=2 wi yip . Since functional F is nonnegative, there exists F ⋆ :=
(ii) There exists a point in R × (−∞, 0) at which inf(a,b)∈P F(a, b). It should be shown that there exists a
functional F defined by (5) attains a value less point (a⋆ , b⋆ ) ∈ P such that F(a⋆ , b⋆ ) = F ⋆ .
p
than ∑n−1
i=1 wi yi .
Before continuing the proof, let us note that Lemma 1.
implies that
Proof. We first prove A(i). Let a(b) := ln( g(ty1 ) ) − bh(t1 ).
1
It is easy to verify that f (t1 ; a(b), b) = y1 and n n−1
F ⋆ < min ∑ i ∑ wi yip .
p

limb→−∞ f (ti ; a(b), b) = 0 for each i = 2, . . . , n. wi y , (6)
i=2 i=1
By definition of the limit there exists a point
b0 ∈ (−∞, 0) such that for all b ∈ (−∞, b0 ), Let (ak , bk ) be a sequence in P, such that
0 < f (ti ; a(b), b) < yi , i = 2, . . . , n. n
∑ wi g(ti )eak +bk h(ti ) − yi
p
F ⋆ = lim F(ak , bk ) = lim

.
Therefore, it follows that k→∞ k→∞ i=1
n n (7)
F(a(b), b) = ∑ wi | f (ti ; a(b), b) − yi| p < ∑ wi yip Without loss of generality, in further consideration we
i=1 i=2 may assume that sequences (ak ) and (bk ) are monotone.
This is possible because the sequence (ak , bk ) has a
for all b ∈ (−∞, b0 ), and therefore claim A(i) holds.
subsequence (alk , blk ), such that all its component
A(ii) The proof is similar to that of part A(i). Let
sequences (alk ) and (blk ) are monotone; and since
a(b) := ln( g(tynn ) ) − bh(tn ). It is easy to check that
limk→∞ F(alk , blk ) = limk→∞ F(ak , bk ) = F ⋆ .
f (tn ; a(b), b) = yn . Since limb→∞ f (ti ; a(b), b) = 0 for Since each monotone sequence of real numbers
each i = 1, . . . , n − 1, there exists a point b0 ∈ (0, ∞) such converges in the extended real number system R, define
that for all b ∈ (b0 , ∞),
a⋆ := lim ak , b⋆ := lim bk .
0 < f (ti ; a(b), b) < yi , i = 1, . . . , n − 1. k→∞ k→∞

Therefore, for every b ∈ (b0 , ∞) we have that Note that −∞ ≤ a⋆ , b⋆ ≤ ∞.


p
F(a(b), b) = ∑ni=1 wi | f (ti ; a(b), b) − yi | p < ∑n−1 i=1 wi yi . To complete the proof it is enough to show that
The rest of the proof is similar, so it is omitted.  (a⋆ , b⋆ ) ∈ P, i.e. that −∞ < a⋆ < ∞ and −∞ < b⋆ < ∞.
Theorem 1. Let P be one of the sets R × (−∞, 0], R × Indeed, the continuity of functional F will then imply that
[0, ∞) and R2 . If the data (wi ,ti , yi ), i = 1, . . . , n, n ≥ 3, are F ⋆ = limk→∞ F(ak , bk ) = F(a⋆ , b⋆ ).

8
8

6 6

4 4

2 2

0.2 0.4 0.6 0.8 1.0 1.2 0.2 0.4 0.6 0.8 1.0 1.2

Fig. 5: Plots of f (t; a(b), b) used in the proof of claim A(i) in Fig. 6: Plots of f (t; a(b), b) used in the proof of claim A(ii) in
Lemma 1. Lemma 1.

c 2015 NSP
Natural Sciences Publishing Cor.
Appl. Math. Inf. Sci. 9, No. 6, 2925-2931 (2015) / www.naturalspublishing.com/Journals.asp 2929

It remains to show that a⋆ and b⋆ are real numbers. To y ✻


20
do this, let us denote
li⋆ := lim (ak + bk h(ti )), i = 1, . . . , n. 15
k→∞

Note that li⋆ 6= ∞ for each i = 1, . . . , n. Indeed, if li⋆ = ∞ 10

for some i, then it would follow from (7) that F ⋆ = ∞,


which is impossible. Also note that li⋆ 6= −∞ for at least 5
one index i because otherwise it would follow from (7)
that F ⋆ = ∑ni=1 wi yip , which contradicts (6). Let index i0 ✲
be such that li0 ∈ R. Then only one of the following three 20 40 60 80 100
t
cases can occur: (i) i0 = 1, (ii) i0 = n, or (iii) 1 < i0 < n.
Now we are going to show that a⋆ and b⋆ are real numbers Fig. 7: The data, initial fit (blue), and optimal fit (red)
in each of these three cases, and so complete the proof of
the theorem. Since li0 = limk→∞ (ak + bk h(ti0 )) ∈ R, note
that it will be enough to show that b⋆ is real. To do this, we
will use the following identities
with F(a0 , b0 ) = 1.966131.
ak +bk h(ti ) = ak +bk h(ti0 )+bk (h(ti )−h(ti0 )), i = 1, . . . , n.
(8) The results of nonlinear LS are given in Table 2
Case (i): i0 = 1. If b⋆ = ∞, then it would follow from (8)
that li⋆ = ∞, i = 2, . . . , n, which is impossible. If b⋆ = −∞,
then it would follow from (8) that li⋆ = −∞, i = 2, . . . , n,
and consequently F ⋆ ≥ ∑ni=2 wi yip , which contradicts (6). Table 2: Least squares parameter estimates and the
Case (ii): i0 = n. If b⋆ = ∞, then it would follow from (8) corresponding value of functional F
that li⋆ = −∞, i = 1, . . . , n − 1, and consequently a⋆ b⋆ F(a⋆ , b⋆ )
F ⋆ ≥ ∑n−1 p ⋆ 27.37 -24.18 1.44361
i=1 wi yi , which contradicts (6). If b = −∞, then
it would follow from (8) that li⋆ = ∞, i = 1, . . . , n − 1,
which is impossible.
Case (iii): 1 < i0 < n. If b⋆ = ∞, then it would follow from The data and graphs of inital and optimal fit are shown
(8) that li⋆0 +1 = ∞, which is impossible. If b⋆ = −∞, then it in Figure 7.
would follow from (8) that li⋆0 −1 = ∞, which is impossible.
Example 2.
Herewith we completed the proof. 
In this example we will illustrate influence of weights
to resulting fit. We will consider the population growth of
5 Numerical examples the United States of America in the period between 2000
Further in the text, we will give few illustrative examples. and 2013. We will use the data from the period between
2000 and 2010 to build an exponential growth model
Example 1. (exponential regression) using few different sets of
weights. The data from 2011 to 2013 will be used to test
Let us consider data set given in the Table 1. how well obtained model fit the data few years in
advance. Let us present the data in Table 3.
This data represent observed average height of
Besides unweighted least squares where all weights are
Cryptomeria Japonica of quality I in the Main Island and
equal 1 (w0 ), we will use following sets of (normalized)
Kiushu (source [24]). To fit this date, we will use
weights:
following two-parameter model:
b
y(t) = ae t . xi
–linear (w1 ): wi = 11
,
Optimal parameters will be obtain in the sense of
(unweighted) least squares, i. e. by minimizing functional
∑ xk
k=1
x2i
29 b –square (w2 ): wi = ,
F(a, b) = ∑ (ae − yi )2 .
ti 11
i=1 ∑ x2k
k=1
As a initial approximation, we use the result of exi /2
linearized model –exponential (we ): wi = 11
,

(a0 , b0 ) = (26.77, −23.46), ∑e xk /2


k=1

c 2015 NSP
Natural Sciences Publishing Cor.
2930 D. Marković, L. Borozan: On Parameter Estimation by Nonlinear Least Squares...

Table 1: The observed average height


Age (ti ) 12 13 14 15 16 17 18 19 20 22 24 25 26 27 28
Height (yi ) 4.1 4.5 4.8 5.4 6.2 6.6 7 8 8.6 9.2 10 10.4 10.6 11.1 11.4
Age (ti ) 30 32 33 35 37 40 46 48 50 56 60 70 75 100
Height (yi ) 12 12.7 13.2 13.5 14 14.8 16 16.5 16.6 17.8 18.5 19.5 19.7 22.1

Table 3: US population (Source: United States Bureau of the Census)


Year (ti ) 2000 2001 2002 2003 2004 2005 2006
Population (yi ) 282162411 284968955 287625193 290107933 292805298 295516599 298379912
Year (ti ) 2007 2008 2009 2010 2011 2012 2013
Population (yi ) 301231207 304093966 306771529 309326295 311582564 313873685 316128839

2
exi Table 4: Least squares parameter estimates and the
–modified exponential (wme ): wi = 11
,
corresponding prediction errors
∑e x2k
a⋆ b⋆ PE
k=1
for i = 1, . . . , 11. w0 -12.8227 0.00923279 7.37169
The results are given in Table 4, where prediction error
(PE) is calculated as w1 -12.8226 0.00923274 7.37054
14
w2
∑ (ea +b ti − yi)2 ,
⋆ ⋆ -12.8225 0.00923268 7.36939
PE =
k=12
we -12.5298 0.00908687 5.48948
for all sets of weights.
wme -10.9354 0.00829343 1.47465
We can conclude that for this particular data, weights
wme give the best prediction. That could be interpreted in
a way that the more recent data have a bigger influence for
the data in the near future. References
The data and graphs of inital and optimal fit with
weights wme are shown in Figure 8. [1] A. Atieg, G.A. Watson, Use of l p norms in fitting curves and
surfaces to data, ANZIAM J. 45(E) C187-C200 (2004)
[2] D.M. Bates, D.G. Watts, Nonlinear Regression Analysis and
Acknowledgement its Applications, Wiley, New York, 1988.
[3] A. Björck, Numerical Methods for Least Squares Problems,
The authors acknowledge the financial support for the SIAM, Philadelphia, PA, 1996.
research by J. J. Strossmayer University of Osijek (project [4] H. E. Burkhart, M. Tomé, Modeling Forest Trees and
“Parameter estimation problem in some two-parameter Stands, Springer, Dordrecht 2012
monotonic mathematical model”). [5] N. Chater, G.D.A. Brown, Scale-invariance as a unifying
psychological principle, Cognition 69 B17-B24 (1999)
[6] C. W. Clark, Mathematical bioeconomics: the optimal
management of renewable resources, Wiley, New York,
y
1990
315
✻ æ
[7] J.E. Dennis, R.B. Schnabel, Numerical Methods for
æ

310
æ Unconstrained Optimization and Nonlinear Equations,
æ
æ
SIAM, Philadelphia, PA, 1996.
305
æ
[8] W.W. Fox, An exponential surplus yield model for
æ optimising exploited fish populations, Transactions of the
300
æ American Fisheries Society 99 80-88 (1970)
295 æ [9] F.R. Giordano, W.P. Fox, S.B. Horton, A First Course in
æ
Mathematical Modeling, Brooks/Cole Publishers, Boston,
290 æ
æ 2014.
æ ✲ [10] Z. Frica, M. Konvicka, Dispersal kernels of butterflies:
æ
2002 2004 2006 2008 2010 2012
t Power-law functions are invariant to marking frequency,
Basic and Applied Ecology 8 377-386 (2007)
Fig. 8: Plots of the data (blue dots), test data (green dots), initial [11] P.E. Gill, W. Murray, M.H. Wright, Practical Optimization,
fit (blue) and optimal fit (red) Academic Press, London, 1981.

c 2015 NSP
Natural Sciences Publishing Cor.
Appl. Math. Inf. Sci. 9, No. 6, 2925-2931 (2015) / www.naturalspublishing.com/Journals.asp 2931

[12] R. Gonin, A.H. Money, Nonlinear Lp-Norm Estimation, Darija Marković


Marcel Dekker, New York, 1989. holds the B.Sc. degree
[13] J. Imaña-Encinas, O.A. Santana, C.R. Imaña, Volumetric in mathematics and computer
and economic optimal rotations for firewood production of science from University
Eucalyptus urophylla in Ipameri, state of Goias, Floresta, of Osijek, as well as
Curitiba 41 905-912 (2011) the M.Sc. in mathematics
[14] N.O. Johnson, A trend line for growth series, J. Am. Stat. from the University of
Assoc. 30 717 (1935) Zagreb. In 2009 she obtained
[15] D. Jukić, A simple proof of the existence of the best a Ph.D. in mathematics from
estimator in a quasilinear regression model, Journal of
the University of Zagreb.
optimization theory and applications 162 293-302 (2014)
Presently she is an assistant professor at the Department
[16] P.A. Marquet et al., Scaling and power-laws in ecological
systems, The Journal of Experimental Biology 208 1749-
of Mathematics, University of Osijek. Her research
1769 (2005) interests cover applied and numerical mathematics,
[17] I. Michailoff, Zahlenmäßiges Verfahren für die Ausführung specifically parameter estimation, mathematical
der Bestandeshöhenkurven, Forstw. Cbl. U. Thar. Fostl. modelling, nonlinear least squares and smoothing
Jarhrb. 6 273-279 (1943) methods. She has co-authored several papers in journals
[18] A.R. Overman, R.V. Scholtz, Mathematical Models of Crop and conference proceedings on these subjects.
Growth and Yield, CRC Press, 2002.
[19] M.J. Panik, Growth Curve Modeling: Theory and
Applications, John Wiley, New Jersey, 2014. Luka Borozan obtained
[20] G.J.S. Ross, Nonlinear Estimation, Springer, New York, his BSc degree at the
1990. J.J. Strossmayer University
[21] F.X. Schumacher, A new growth curve and its application to of Osijek, Department
timber-yield studies. Journal of Forestry 37 819-820 (1939). of Mathematics. Currently
[22] G.A.F. Seber, C.J. Wild, Nonlinear Regression, John Wiley, he is a MSc student at
New York, 1989. the same university. His
[23] A. Spirling, The next big thing: research interests are
scale invariance in political science, numerical computation
http://www.people.fas.harvard.edu/∼spirling/documents/powerlawSend.pdf
and software engineering.
[24] W. Terazaki, Notes on the Analytical Interpretation of
Growth Curves for Single Tree and Stands, and on
Application for the Construction of Yield Table for
Cryptomeria Japonica, Bulletin of the Forest Experiment
Station, 151-202 (1915)
[25] D. Zajdenweber, Scale Invariance in Economics and in
Finance, Scale Invariance and Beyond 7 185-194 (1997)

c 2015 NSP
Natural Sciences Publishing Cor.

You might also like