Abstract
In many mathematical models with complex parameter functions, it is critical to suggest a strategy for estimating those functions. This paper considers the problem of fitting parameter functions to given data by using a finite number of piecewise functions. Here, the primary aim is to propose several necessary and sufficient conditions that can guarantee the existence of parameter functions. With these criteria, the paper provides an algorithm for solving the parameter estimation problem efficiently. Numerical examples verify that the proposed algorithm is effective in estimating parameter functions in differential equations.
1 Introduction
Mathematical modelling using a system of differential equations with a set of parameters has been widely used in most scientific fields, including physics, chemistry, biology, earth and environmental sciences, and engineering. Mathematical problems related with such parameter-dependent differential equations can be classified into forward and inverse problems. Forward problems need the analysis and simulation of the model for given parameter values, whereas inverse problems require the estimation of parameters based on the measurement of output variables. The latter have received considerable attention in the last three decades (see e.g. [Citation1, Citation2] and the references therein). As statistical approaches to parameter estimation problems, nonlinear regression and other data-fitting methods have been widely used.[Citation3–Citation7]
Among various important problems related to parameter estimation, this paper is particularly interested in the problem of fitting parameter functions to given data by using a finite number of piecewise functions. The primary aim of the paper is to propose some conditions for the existence of parameter functions and introduce an efficient algorithm for estimating these functions. The number of subintervals and the type of basis function in each subinterval are not known a priori, and therefore the main concern is how the number of subintervals and the basis function can be determined and how time nodes can be optimally found for each subinterval.
This paper begins with linear stage-structured differential equations whose parameter functions can be approximated by piecewise constants, and other types of basis functions for approximating parameter functions and nonlinear differential equations with noisy data are left to future research. For a given data-set, the focus is on finding a piecewise constant parameter function in a differential equation to interpolate the data. A necessary and sufficient condition for the existence and uniqueness of the piecewise constant parameter function when subintervals are known is addressed first. In addition, several criteria for determining time nodes and then parameter functions when subintervals are yet to be determined are proposed. These criteria are provided by a -test function in Section 2. An algorithm is suggested to determine time nodes and coefficients for each subinterval optimally with the given data-set and the number of subintervals. Several numerical experiments demonstrate that the algorithm is effective in estimating parameter functions. Further, the solution of the linear stage-structured differential equation with the estimated parameter function provides a good approximation for the solution of the nonlinear differential equation.
The rest of this paper proceeds as follows: Section 2 starts with a formulation of the problem of estimating parameter functions with a linear stage-structured differential equation. A -test function is defined to determine whether there exist suitable solutions of the differential equation to interpolate a set of given data. Section 3 presents the necessary and sufficient conditions for the given data-set to solve the problem of estimating parameter functions when time nodes of subintervals are known. Section 4 derives the conditions by using the -test function for determining parameter functions when time nodes of subintervals are yet to be determined. Section 5 introduces an algorithm for determining parameter functions based on an optimization problem, and Section 6 presents the numerical results for estimated parameter functions and their corresponding solutions based on the proposed algorithm.
2 Parameter estimation problem
2.1 Motivation
Numerous biological phenomena can be described by stage-structured modelling. For instance, growth rates of most fish species are sensitive to seasonal effects, mainly to temperature changes, and it is well known that there is a strong relationship between specific growth rates and temperature for most populations of bacteria (see e.g. [Citation8–Citation10] and the references therein). In addition, many delay differential equations with a stage structure can be found in the literature.[Citation11–Citation15]
For the given data-set the following stage-structured differential equation is considered:with the parameter function1 1 where the number of intervals , the node and the parameter , are to be determined by assuming that Here, , and depend on the differential equation and the given data-set .
2.2 Linear case
This paper focuses on linear stage-structured differential equations by emphasizing the relationship between the parameter function and the data-set. Nonlinear differential equations are left to future research.
Now consider the following estimation problem: For a set of given data find in the form of (Equation11 1 ) such that2a 2a 2b 2b
Note that the solution of (2) can be given by3 3 when the parameter function in (Equation11 1 ) is known.
Because most physical data have positive values, solutions of Equation (Equation2a2a 2a ) are assumed to be positive and data-sets are restricted to cases such as those described in Definition 2.1.
Definition 2.1
The data-set is called positively well-posed if has the following two properties:
(i) | is a strictly increasing sequence, | ||||
(ii) |
If there is only a single subinterval, then the necessary and sufficient condition for the constant function to be uniquely determined in the interval is summarized in the following lemma:
Lemma 2.2
Suppose that there is a given data set that is positively well-posed in a single interval . Then the constant function and the solution satisfying (2) can be uniquely determined if and only if
Proof
Suppose that satisfiesThen the given data-set must satisfyThat is, are collinear. The existence and uniqueness of and are equivalent to
Before describing the relationship between data points, the -test function is first defined.
Definition 2.3
Define the -test function by
Proposition 2.4
Suppose that . If the -test function C is positive (or negative), then there exists a log-concave (or log-convex) function that passes through , and . In addition, the following properties hold:4 4
Proof
Note that the -test function can be rewritten asSuppose that . Then which implies the existence of a concave quadratic function that passes through and The case of is similar. Equation (Equation44 4 ) follows immediately from the definition of the -test function .
Remark 2.5
implies that and lie on one exponential curve: whereEquivalently, the exponential curve can be rewritten as
Definition 2.6
The data-set is called reducible for the parameter estimation problem if and only if there exist three adjacent data points such that Otherwise, the set of data is called irreducible.
Remark 2.7
Suppose that the given data-set is positively well-posed and irreducible in the interval with the partition . Then every subinterval contains no more than two data points if the data-set determines the parameter function to be uniquely fulfilling (2).
3 The case of known subintervals
This section starts by looking at the case in which (), are given. The following theorem provides the conditions for the given data-set to guarantee the existence and uniqueness of the parameter function and the corresponding solution of (2).
Theorem 3.1
Let Suppose that the given data set is positively well-posed and let be given time nodes of Then there exists a unique parameter function and its solution that satisfy (2) if and only if the following condition holds:
If any interval contains two data points, then data points are located in for any .
Proof
Rearrange data points as follows:
(i) | For , let be the number of data points in the interval so that and , | ||||
(ii) | For , denote by the set of data points in the interval , | ||||
(iii) | For , is a strictly increasing sequence. |
With some elementary operations, can be transformed into the following matrixwhere denotes a zero block matrix andandSincethe following can be obtained by recalling :Based on the above, in the case of , , and thus is not invertible. Therefore, there is a need for the condition with to make invertible.
Remark 3.2
(i) | Theorem 3.1 provides the necessary and sufficient condition for data and subintervals for the determination of a parameter function uniquely to fulfil (2). That is, if a subinterval contains two data points, then the number of data points should equal that of subintervals in the whole left side of the interval. | ||||
(ii) | In the case of and therefore is not invertible. Although the parameter function can be determined, the uniqueness of and its corresponding function satisfying (2) cannot be guaranteed. | ||||
(iii) | In the case of there is an over-determined problem. The only case in which the over-determined system has a unique solution is that the linear system can be reduced to the system with |
Remark 3.3
Let a given data-set that is positively well posed and be given time nodes in with Suppose that there are two consecutive subintervals and containing two data points in each subinterval. Then the parameter function cannot be determined uniquely. Indeed, for to be determined uniquely, Theorem 3.1 requires that and contain and data points, respectively, which leads to a contradiction.
4 The case of subintervals not yet determined
The previous section describes the necessary and sufficient condition for a data-set to determine a parameter function and its corresponding solution uniquely when The parameter function estimation problem is over-determined in the case of If the problem is over-determined with too much information, then a particular solution may not exist. As one way to address this problem, the number of unknown parameters can be increased by assuming that time nodes need to be determined. Suppose a given data-set is positively well-posed and irreducible. Then every subinterval should have at most two data points. Note that in each subinterval, the constant value of the parameter and its corresponding solution is uniquely determined by two data points. Therefore, to determine the parameter function for the whole interval, the following types of contiguous subintervals need to be considered:
(T1) | The first and last subintervals in this category both contain two data points, whereas all other subintervals contain only one data point, | ||||
(T2) | Only the first or last subinterval in this category contains two data points, whereas all other subintervals contain only one data point. |
Lemma 4.1
Assume that and are any nonzero real numbers. Then the following two conditions are equivalent:
(i) | or , | ||||
(ii) | or |
Proof
Suppose (ii) is not true, that is, Then which implies that (i) is not true. By contrapositive, this completes the proof.
Lemma 4.2
Suppose that and the given data set is positively well-posed and irreducible over the whole interval Then there exist unique and the corresponding parameters and such that the function given by (Equation33 3 ) satisfies (2) if and only if
Proof
Suppose that there exist and the corresponding parameters and such that the function given by (Equation33 3 ) satisfies (2). Then we haveand (Equation2b2b 2b ) also requires the following relationships:The parameter values of and and the logarithmic values of and can be determined in terms of by solving the above equations:7 7 From the continuity condition , we have8 8 Combining (Equation77 7 ) and (Equation88 8 ) gives9 9 and finally, Now the condition in which lies between and needs to be derived (see Figure ). This condition can be restated as and a direct calculation givesNote that the denominator is not zero because data points are irreducible. Therefore, if and only if
Lemma 4.3
Suppose that and the given data set is positively well-posed and irreducible over the whole interval Then there exist and corresponding parameters and such that the function given by (Equation33 3 ) satisfies (2) if and only if
Proof
Suppose that there exist two free nodes and such that following conditions hold (see Figure ):10 10 For this problem, we need to determine eight variables and . From (Equation1010 10 ), it is clear that11 11 By the continuity condition of the solution12 12 Therefore, and can be determined by and the data-set . However, is related to and , and therefore showing the existence of is equivalent to showing that the algebraic systemhas at least one solution . This leads to the relationshipwhich is equivalent to13 13 Combining Equations (Equation1111 11 ), (Equation1212 12 ) and (Equation1313 13 ) gives14 14 The RHS of (Equation1414 14 ) is denoted as . Then the roots of the function in the rectangular domain are considered. Since is a polynomial of and the function has roots in the domain provided thatRepresenting the above conditions by using the -function (see the Appendix) giveswhich completes the proof.
Lemma 4.3 can be extended to general cases.
Theorem 4.4
(-test) Suppose that ( is a positive integer) and the given data set is positively well-posed and irreducible over the whole interval Then there exist for all , and corresponding parameters such that the function given by (Equation33 3 ) satisfies (2) if and only if15 15 for some
Proof
As in Lemma 4.3, the theorem holds for . Assume that it also holds for .
()Suppose that there exist for all , and corresponding parameters such that the function given by (Equation33 3 ) satisfies (2). By the continuity condition as in Equation (Equation99 9 ), we haveNow consider the new data-set Since the theorem holds for the following condition holds:16 16 for some . Define by that is,Since is a linear function of and , the condition (Equation1616 16 ) impliesFrom the calculation, we getandTherefore,By Lemma 4.1, we obtainThis means that there exists such that (Equation1515 15 ) holds.
()Suppose that (Equation1515 15 ) holds for some Then consider the new data-set Since the theorem holds for there exist for all , and corresponding parameters such that the function given by (Equation33 3 ) satisfies (2). If is chosen, then can be determined such that the function in the subinterval passes the data point Suppose that (Equation1515 15 ) holds for some Then consider the new data-set By similar arguments, can be constructed for all , and corresponding parameters
Remark 4.5
Under the assumption of the irreducibility of data points, the value of the -function is always nonzero. From this fact, the condition (Equation1515 15 ) can be rewritten as (Equation1717 17 ).
Theorem 4.6
(modified -test) Suppose that ( is a positive integer) and the given data set is positively well-posed and irreducible over the whole interval Then there exist for all , and corresponding parameters such that the function given by (Equation33 3 ) satisfies (2) if and only if17 17 for some
Then the existence of a solution of Equation (2) can be easily checked by the above theorem when . However, this theorem cannot guarantee the uniqueness of , the parameter and the solution In this regard, a meaningful algorithm is now introduced to determine the time nodes, the parameter function and its corresponding solution by formulating an appropriate optimization problem.
5 An algorithm for estimating parameter functions
This section introduces an algorithm for estimating a piecewise constant parameter function in the case of If input information (data points, the number of subintervals and the structure of subintervals) is given, then the whole interval can be divided into several subgroups that are of the form T1 or T2. As discussed earlier, the modified -test is conducted only on T1 because the existence of and is always guaranteed in the case of T2. If the modified -test fails, then the existence of and cannot be guaranteed, and therefore another structure of subintervals needs to be chosen.
Assume that the modified -test is successful in T1. Then can be found by minimizing the total variation of the constant parameter in the whole interval. For this, the objective functional can be defined as18 18 where is a constant parameter in the subinterval depending on Here, is the number of intervals in T1 or T2. In addition, there are following inequality constraints for the form of T1(see Figure ):19 19 Since can be expressed uniquely in terms of and the domain variables of the functional (Equation1818 18 ) can be reduced by one. Then the explicit formulation for the functional can be rewritten as follows:whereandfor . Combining these recursive relationships givesfor This implies that the term can be represented by time nodes and data points.
In a similar manner, the constraint in the condition (Equation1919 19 ) can be expressed for the variables asThis is obtained by applying Lemma 4.2 to the following set:Here, the numerical method proposed in [Citation16] is employed to find that can minimize the objective functional subject to inequality constraints. Now the proposed algorithm is briefly introduced (see Algorithm 1).
6 Numerical simulation
6.1 Data-set 1
Data-set 1 deals with a set of 10 data points listed in Table .
These data points are generated from the analytic solution to the following logistic differential equation.20 20 where the parameters , and are from [Citation17]. Note that this is a positively well-posed data-set.
Here, the goal is to find a piecewise parameter function for the stage-structured differential Equation (2) that minimizes the functional (Equation1818 18 ). Before proceeding, choose an appropriate number of subintervals and the structure of the subinterval. Because 10 data points are considered, the appropriate number of subintervals, , ranging from 5 to 9. If is less than 5, then it goes against to the definition of an irreducible data-set. If is greater than 9, then its subinterval contains only one or zero data point, and this does not properly reflect on the dependency of data points. For this reason, select .
In this example, assume that the number of subintervals is seven, that is, , and the structure of subintervals is (2, 2, 1, 1, 1, 1, 2), which means that the first, second and seventh intervals have two data points, whereas the others, one data point. As shown in Figure , the whole interval is composed of two subgroups (G1 and G2) that are rectangled to identify which data points belong to which subgroups. This decomposition is based on Lemma 2.2 and determines parameters uniquely by using two data points in a single subinterval.
First, consider the subgroup G1. Note that this subgroup has two subintervals and is classified as T1. According to Algorithm 1, the modified -test needs to be conducted for G1. From the simple calculationthere exist and for the solution of Equation (2), as shown in Lemma 4.2.
Similarly, the second group G2 has six subintervals and is classified as T1. Here, of course, two data points and are common elements of G1 and G2. The subgroup G2 is also T1, and therefore the modified -test is again needed. Sincethe existence of time nodes and parameters in G2 is guaranteed. Then the results shown in Table can be obtained by applying a minimization scheme to minimize the functional (Equation1818 18 ).
With estimated time nodes and parameters, the error between the derived and exact solutions of (2) is . Figure shows the graph of the solutions.
6.2 Data-set 2
Data-set 2 generates as random numbers in and by solving the following differential equation at :21 21 The generated data-set is shown in Table . Data-set 2 considers seven subintervals (), and the structure of subintervals is (1, 1, 1, 2, 2, 1, 2). In other words, the fourth, fifth and seventh subintervals have two data points, whereas the others, one data point. As depicted in Figure , the whole interval can be decomposed by three subgroups G1, G2 and G3. The first subgroup G1 has four subintervals and is classified as T2. In this case, there is no need to conduct the modified -test because it is T2, and therefore the existence of time nodes and corresponding parameters are always guaranteed. On the other hand, G2 and G3 have two and three subintervals, respectively, and both are classified as T1. Therefore, the modified -test needs to be conducted. In G2, however, the value ofis negative, implying the nonexistence of . Accordingly, it is not possible to determine a piecewise constant parameter function satisfying (2) with this structure.
Consider another structure of subintervals. Suppose that seven subintervals () are given and the structure of subintervals is (1, 1, 2, 1, 2, 1, 2). In this case, the third, fifth and seventh subintervals have two data points, whereas the others, one data point. The whole interval can be decomposed into three subgroups G1, G2 and G3. As in the cases of the other data-sets, the modified -test needs to be conducted for G2 and G3, which are of T1 type. Here simple calculations indicate that the modified -test is successful, and therefore the optimal parameter function can be calculated. Table shows the results. In addition, the error between the derived and exact solutions of (Equation2121 21 ) is Figure shows the graph of the solutions with estimated time nodes and parameters.
7 Funding
This research was supported by NRF of Korea [grant number 2011-0000344], [grant number 2009-0065241], [grant number 2008-C00043].
References
- Beck JV, Arnold KJ. Parameter estimation in engineering and science. New York (NY): John Wiley & Sons; 1977.
- Tarantola A. Inverse problem theory. Methods for data fitting and model parameter estimation. Amsterdam: Elsevier; 1987.
- Banks HT, Bihari KL. Modelling and estimating uncertainty in parameter estimation. Inverse Probl. 2001;17:95–111.
- Banks HT, Fitzpatrick BG. Statistical methods for model comparison in parameter estimation problems for distributed systems. J. Math. Biol. 1990;28:501–527.
- Hvala N, Zec M, Strmc̆nik S. Non-linear model parameter estimation-estimating a feasible parameter set with respect to model use. Math. Comput. Modell. Dyn. Syst. Meth. Tools Appl. Eng. Related Sci. 2008;14:587–605.
- Martinsons CD. Optimal parameterization of a mathematical model for solving parameter estimation problems. Inverse Probl. Sci. Eng. 2005;13:109–131.
- Seber GAF, Wild CJ. Nonlinear regression. Hoboken (NJ): Wiley-Interscience; 2003.
- Beverton RJH, Holt SJ. A review of the lifespans and mortality rates of fish in nature, and their relation to growth and other physiological characteristics. Chichester: John Wiley & Sons Ltd; 2008. p. 142–180.
- Cren EDL. The length-weight relationship and seasonal cycle in gonad weight and condition in the perch (perca fluviatilis). J. Animal Ecol. 1951;20:201–219.
- White P, Kalff J, Rasmussen J, Gasol J. The effect of temperature and algal biomass on bacterial production and specific growth rate in freshwater and marine habitats. Microb. Ecol. 1991;21:99–118.
- Aiello WG, Freedman H. A time-delay model of single-species growth with stage structure. Math. Biosci. 1990;101:139–153.
- Aiello WG, Freedman HI, Wu J. Analysis of a model representing stage-structured population growth with state-dependent time delay. SIAM J. Appl. Math. 1992;52:855–869.
- Gourley SA, Kuang. Y. A stage structured predator-prey model and its dependence on maturation delay and death rate. J. Math. Biol. 2004;49:188–200.
- Gurney WSC, Nisbet RM. Fluctuation periodicity, generation separation, and the expression of larval competition. Theoret. Population Biol. 1985;28:150–180.
- Nisbet RM, Blythe SP, Gurney WSC, Metz JAJ. Stage-structure models of populations with distinct growth and development processes. IMA J. Math. Appl. Med. Biol. 1985;2:57–68.
- Rockafellar RT. The multiplier method of Hestenes and Powell applied to convex programming. J. Optim. Theory Appl. 1973;12:555–562.
- Gauze GF. The struggle for existence. Baltimore (MD): The Williams & Wilkins company; 1934.
Details of calculation in Lemma 4.3