Kárný, Miroslav

English

Institute of Mathematics AS CR, v. v. i.

KybernetikaMiroslav KárnýRecursive parameter estimation of regression model when the interval of possiblevalues is givenKybernetika, Vol. 18 (1982), No. 1, 37--49Persistent URL: http://dml.cz/dmlcz/125590Terms of use:© Institute of Information Theory and Automation AS CR, 1982Institute of Mathematics of the Academy of Sciences of the Czech Republic provides access to digitizeddocuments strictly for personal use. Each copy of any part of this document must contain theseTerms of use.This paper has been digitized, optimized for electronic delivery and stamped withdigital signature within the project DML-CZ: The Czech Digital Mathematics Libraryhttp://project.dml.czK Y B E R N E T I K A - V O L U M E 18 (1982) , N U M B E R 1 RECURSIVE PARAMETER ESTIMATION OF REGRESSION MODEL WHEN THE INTERVAL OF POSSIBLE VALUES IS GIVEN MIROSLAV KARNY The paper presents aposteriori maximum likelihood estimation of coefficients of regression model when prior bounds on their possible values are given. Its objective is twofold: — conceptually clarify the problem in which the difference between experience accumulation and decision-making became apparent — to describe modification of Fletcher-Jackson algorithm [1] which admits recursive parameter estimation in apriori given interval. 1. INTRODUCTION The main feature of the Bayesian approach to the system identification [2] is its capability to feed prior information into the problem. The simplest but important prior information is that about the range S of possible values of the unknown parameters 0. This type of information can be, in principle, introduced in a very simple way: The support of the prior probability density function (which exists in the case assumed) p(&) coincides with S i.e. (1) p(0)<x<?(0,O)Xs(0) where the symbol oc means proportionality, the function x s(0) is the indicator of the set S and i£(0, 0) is a nonnegative function for which (2) 0 < / ( 0 ) = Җ , 0) d < oo . It is known, however, that only relatively narrow range of identification problems can be exactly converted into practically feasible algorithms. The one of important exceptions is the linear-in-parameters multivariable normal regression model (briefly fundamental model) [3] with conjugate prior distribution 37 function [4], [5]. In this case the Bayesian identification leads to the algorithm which formally coincides with the (recursive) least squares. This coincidence often results into the opinion that Bayesian approach is something superfluous, at least for the set of linear-in-parameters models. But assuming the case of the fundamental model with bounded range of possible values of parameters we arrive to different solutions using different interpretations of the standard least squares (different projection see e.g. [6], [7]). We shall try to clarify the source of this diversity. It can be easily found that in the bounded range case also for fundamental model the full Bayesian solution is algorithmically unfeasible. There are at least two ways to overcome this trouble: To select conjugate prior distribution which is practically equal to zero out of the range S. This is often satisfying solution. But as an apriori unknown value of the likelihood function is decisive for what can be taken as zero probability event, this solution may be dangerous, when point estimates are needed. To restrict oneself to the case of point estimates maximizing the aposteriori likeli-hood function. The main part of this paper deals with this problem. 2. PROBLEM FORMULATION AND PRELIMINARIES Let us assume that a controlled system can be described by multivariate linear-in-parameters normal regression model (briefly fundamental model) (j) p(y(t)\S"l\u(t),e) = N(y(l),Q'1) where we denote te {1, 2 , . . . } — discrete time p(a | 0) - probability density function of a. conditioned on /5 y(t) - v-vector of outputs u(t) — ^-vector of inputs d(t) = (y(t), u(l)) - measured data dw = (d(1), d(2), •-., d(t)) - process history N(y, Q-1) = ~ 1/2 exp {-\(y - yf Q(y - y)} 2K y(t) = PTz(t) — regression function z(t) = z(d(t~1), M(o) ~ known g-dimensional function © = (P, Q) — unknown system parameters 38 We assume that (2) & e S = {(P, Q) : P g, P jg P, Q > 0} where the given boimds P, P satisfy (3) - o o < P < P < + o o and vector inequalities in (2) and (3) are interpreted entrywise. The inequality Q > 0 means that Q is a positive definite matrix. Possible controllers are restricted to the class which fulfils natural conditions of control [2] i.e. explicit values of unknown parameters cannot be used K«cok'-^0) = K"(oK~ir)-Assuming conjugate prior distribution restricted to S, we obtain the following solution of the Bayesian system identification: the aposteriori probability density function of unknown parameters is given by (4) P( \d") = *шмň I(0 where ("conditional") likelihood function evolves according to (5) <£(0, t) = 3>{e, t - 1) p(y(t) [ d«"1} , u(l), 0) ^(0, 0) is given in (IT), (1.2) and (6) / ( 0 = (&(e, t) Xs(0) d o = f <£(0, t) d0 Bayesian prediction is given by (7) p(y(t)\d^\u(t))^^-Ht-i) Specializing relations (4) to (7) for the system model (1) we have W-M4~/H~/])}= exp {-I tr [Q(PTVm P -PTKm - Vjy(l)P + V„(0)]] = 2 exp { - I tr [Q((P - P(t))T Vm (P - P(t)) + A(t))]} 9(0 = »( , -„ + 1 (8) ^ ( , í) = Q 2ҡ Q 2łt = Q 2л where (9) 39 (10) V(0 -о = \VM KW!) v = V D + M M |Л*0 V«) J}_ UoJUoJ (11) í w = S (12) é(0 = .V(t) - PJt-i^t) (13) A(I) = Fy(0 - VTmVz-(t)Vzm = A((_1} + 1 + z(t)C(t^г)z, (0 / 1 4 \ C _ y-i _ C _ e(.-l)Z(Q2(0e(.-l) V1 J (0 - Vz(t) - L (r- l) T i + -(,)-•((_ i)Z(t) Recursive parts of the relations (11), (13), (14) are implied by the relation (10). The entities (9) and (10) form sufficient statistic of the problem assumed, their recursion results from Eq. (5) and initial conditions from Eq. (1.1). The condition (1.2) is satisfied for S ( 0 ) _ 0, V(0) > 0. The form of the recursive least squares (11), (12), (13), (14) is nothing more than the alternative recursion for equivalent sufficient statistic (9, P, A, C). It implies immediately the following conclusions: . Information contained in measured data is compressed in statistic (9, P, A, C) and any prior information cannot destroy (i.e. singularity transform) this statistic. . When the statistic PU) is the admissible estimate of P i.e. (15) P _ P(t) ^ P then maximum likelihood estimate PML(t) coincides with this "least square" estimate. . The prediction (7) is always influenced by bounds P, P, even when (15) is fulfilled. Apparently computation of the predictive probabihty density (7) is a difficult problem which can be practically solved only for small g [8] and will be omitted here. . Maximum likelihood estimates of the parameters 0 are given: (16) QML = 9[(PML - P)T VZ(PML -P) + A]"1 (17) PML = arg min (tr [QML(PTVZP - PTVzy - VTP)]} P&P^P It follows that for multi-output systems (v > l) the estimation of Q and P cannot be separated as in unbounded case. Probably some sort of bootstrap construction can be successfully used, but we restrict ourselves to the single-output case (v = l) in which such splitting appears, i.e. (18) PML - arg min [i(PTVzP - 2PTVzy)] . pgpgp Algorithmic solution of the problem (18) is contained in the next part. 40 3. FLETCHER-JACKSON ALGORITHM Searching for maximum likelihood estimates (v = 1) we solve the sequence of the following problems (1) min [i(P rVz( t )P - 2V^(r)P)] r - 1,2,... We assume that the regular prior distribution is selected and sufficient condition for validity of (1.2) is fulfilled (2) V2(o)>0, 3 ( 0 ) >=0. It follows that (3) Vm > 0 for all t. Under the assumption (4) - oo < P < P < + oo the next theorem can be easily proved. Theorem 1. The solution of the problem (l), (3), (4) exists and it is unique. Let us describe briefly the algorithm, which Fletcher and Jackson [ l ] propose for solution of problem (l), (3), (4) in the fixed time r. The theoretical behavior of this algorithm [ l ] comprimes Theorem 2. The algorithm described below needs finite number of iterations. Algorithm. 1) Set r = 0, select initial admissible approximation (5) P g P r <\ P and compute gradient of the minimized function (6) g(F) = V2F - Vzy . Sort P r in two parts (7) F = ppl} k Q where (8) P; < P,. < P; and g(F) = 0 for i = 1, 2 , . . . , k (free variables) (9) P,. = P,. or P; = P, or P, < Pi < P; and g{p) + 0 for i = fe + 1, fe + 2 , . . . , Q (fixed variables) 41 2) Select some fixed variable Ph 1 > k which will be changed by a step a. Admissible signs of changes are apparently restricted (10) a > 0 if P ; = P, a < 0 if P ; = P ; We take such step size that free variables (8) can be changed in admissible set in such a way that (11) g{Pr+1) = 0 for i = 1,2, . . . ,* It can be easily found that we can achieve a smaller value of the minimized function iff for admissible a (7,) (8) such / exists that (12) gt(F) a < 0 i.e. if such / and admissible a does not exist optimal solution is achieved. The optimality conditions described above are apparently specialized form of Kuhn-Tucker conditions. 3) If we find some / which does not fulfil Kuhn-Tucker conditions we change P r to F+1 according to the rules (13) P r + 1 = PI + a (14) P r + 1 = P r i> k, i 4= / (15) P- + 1 = P- + ayi i = 1, 2, ..., k where /(-vector y is determined from validity of the Eq. (11): The step size a must result in an admissible P r + 1 . There are the following possible branches: i) The step a0, resulting from requirement a ;(Pr + 1) = 0 is admissible. Its selec-tion moves the entry Ft + 1 into the part 1P, the dimension * is increased, ii) The a0 is inadmissible, but there is admissible step which transfer l-th entry to the opposite boundary and 1pr+1 does not cross admissible bound. This step size does not change dimension of *P but the entry P r + 1 now apparently fulfils Kuhn-Tucker conditions, hi) Also the step size from ii) is too long. Some free variable, say Pm+1 , m is k, achieves its bound. Using this step size we start "minor" iterations. We do not change the selected index / and try to improve the approximation changing this entry. Only the variable Pm is moved from the part iP to the 2P, the dimension k decreases. These minor iterations must finish at most in k ^ Q steps because always the entries 1P are moved to the opposite part 2P. 4) Re-computing P r to Pr+l we re-compute also g(Pr) to o(P, + 1) and return to the point 2); only in minor iterations (the case iii) of the 3)) we return to the point. 3). 42 For effective and numerically stable computation of the vector y (15) Fletcher and Jackson propose to use partial LDLT factorisation. They assume V, decomposed (16) V = LDLT =- r iv i2Vr]} k = r IL o i r 1 ^ o i p u I 2LTT k ii2v 2VJ L 1 2 L 2 L J L ° 2 D J L ° 2Lr_ T /T where Lis lower triangular matrix with units on the main diagonal, D is diagonal matrix. The partition in submatrices corresponds to the partition (7). Partial decomposition means that the submatrices lL, 12Land 1D can be computed without computation of the 2L, 2D. Moreover the following identity is used (17) 2H = 2L2D2LT = 2V-l2V'Dl2VT With notations (16) and (17) we can find that the vector (15) solves triangular system (18) >LTy = - c o l , ^ 1 2 ^ ) and the following recursion is valid for nonzero part of the gradient (19) 9i(F+1) = 9i(F) + a 2 / / , . -M -„ i = k + 1 , k + 2,...,g Partial decomposition can save computational effort, but in the problem discussed it has two disadvantages . the partial decomposition of V(() cannot be efficiently used in the sequence of problems i.e. in transformation of minimized function by new data. . the computation in Eq. (17) may be numerically unstable: positive definite matrix 2H must be obtained as a result of substraction. These two points forced us to use full LDLT decomposition. In the sequence of problems it does not increase computation effort substantially and it removes pos-sible numerical troubles. Now we describe in which way this modified Fletcher-Jackson algorithm can be implemented into data accumulation. It seems to be reasonable to take as an initial approximation at time t the optimal estimate obtained in the time t — 1 i.e. (20) IML(r) = 'ML(t-l) The relations (2.10), (2.12) and the choice (20) result in the following recursion for the initial gradient: (21) 9(,)(P°ML(,)) = 9(t)(PML(t-i)) = VoiWi) - Ky(» = = Vl(l-l)PML(t-l) - Ky(t-l) ~ Z(t)(y(t) - FML(t-l)Z(t)) = = 9(t-l)(PML(t-l)) - Z(tf(t) 43 Having LDLT decomposition of the Vz(,_i) the recursion (10) for the V, can be written ( 2 2 ) ht)D(t)L\t) = L(t-i)(D(t-i) + f(t-i)fT(t-i))Llt-i) where the vector / (superfluous index t — 1 will be suppressed) solve the triangular system (23) Lf = z ( t ) Using definition of the LDLT decomposition we decompose (24) D + ffT = QD(t)QT where the Q is the matrix of the same type as the L. Then (25) _«,,=-v..,g Algorithm which solve (23), (24), (25) can be written in the following compact form (vector z ( ( ) is destroyed) (26) S = 1 j = l , 2 . . . . , в C = a = b = z y G,v 1/s s = s + c * * 2/ö a = ö * s GJJ = a * Ь a = cja b = j + 1J+2,.. GІJ ., g 2 | = zt — b * c G,v = b + a * 2, end of i, / It means that the use of the full LDL1 decomposition requires approximately only two times more multiplications than the simple re-computation of Vz according to (2.10). It requires aproximatelly Q2 multiplications. The last problem which has to be solved is reconstruction of the full LDLf decom­position after permutations performed in the steps 3i), 3iii) of the Algorithm. As this problem seems to us to be of independent interest (the task appears in some classification problems), the solution will be given separately. 44 4. THE PERMUTATION OF VARIABLES AND LDLT DECOMPOSITION We have given LDLT decomposition and the quadratic form (1) ijj(x) = xTLDL'x The vector x splits in two subvectors (2) m We solve the following problems: a) To remove the variable x,„, m <_ k, from lx. For this purpose we exchange (3) *m*-* xm+\ then we transform the destroyed LDL1 decomposition into right form, symbolic-ally we perform transformation £.. We repeat procedure up to /c-th entry, i.e. (4) xm *->xm+10l x,„ <-> xm + 2 0$ ... xm <^>x_0l then we set k = k — 1. b) To transfer /-th variable, / > k, from 2x to 'x. With notations of the problem a) we proceed (5) x, <-» x,_i £#x,<-» x,_2 5? ... x, <-> x t + i 5? and then we increase fe to k + 1. The solutions of the both problems require effective computation (6) x,<->x i + 1_? Exchanging the i-th and (i + l)-st variable we exchange the i-th and (i + l)-st rows in the matrix L. This permutation results in the following forms of the i-th and (i + l)-st columns (7) - 1 { "o " 1{ a 1{ 1 л . - 1 { 1{ 0 - i + 1 1{ e - / + 1 The new i-th and (i + l)-st columns must have the form t -1{ "o" 1{ 1 1{ b J_ i-í{ 1{ . - / + 1 1{ - i + 1 The original i-th and (i + l)-st terms of the diagonal matrix D, say _i, _2, we transform in rx, r2 45 The following identity must be valid (9) a 1 " 1 0 Лh 2Һ Гqi0Ҡa 1 t й П - ф O І Г r , O І Г l Ь / " [0 q2Jll 0 2/2тJ Ь 1 [0 r 2 J [ o 1 2 /т l j 2 j The equation (9) is solved by (ю) r. = a2qi + q2 Ь = ^ -Гl r = îlîl Гí J = yhb + 2Һ 2 j = III - 2lîЯ <72 i.e. the transformation (6) requires aproximately 3(g — i — l) multiplications and the J 2 - 1 sequence of the type (4) or (5) between variables j , < J 2 needs 3 __ (g — i — 1) i-Ji multiplications. This quantity is bounded by 3/2g2 multiplications. 5. ILLUSTRATIVE EXAMPLE To illustrate the behavior of the resulting algorithm we give rather simple example of parameter identification of a controlled system in the open loop. The system is determined by (D ý{() = p" Z ( f ) = [2, 0.98, -1 .98] Уit-n , 0 = 1 and the white-noise input is normally distributed (2) « ( 0 cr JV(0, 1) The estimates of P 2 = 0.98 are drawn in the Fig. 1 for unbounded and bounded case in which (3) P = 0" , P = T 0 1 - 2 0 Because the true parameter values are in the range given both estimates coincide "asymptotically". But transient identification part is influenced substantially as is apparent from the sample dispersion of prediction errors which is 5.6 for unbounded and 1.2 for bounded case. 1.51 1.0 I I.5І 4.21 J ^ ^ ^ V v / ^ ^ ^ jnbounded case 10 20 30 40 50 t ғ i g - L The feature that the activity of any bound influences all estimates is illustrated in Fig. 2 which contains the estimates of P2 when bounds on P 3 were changed to the unrealistic ones (Pt « 2.1, P3 = 0) (4) 0 ^ P 3 ^ 2 1.0-0 8 -O.6.. 0.4.. O.2.. 0 Estimates of r| when bounds on E>ML are unrealistic ( ^ « 2.1 , §ni_=0) 47 6. CONCLUSIONS The Bayesian solution of the problem of the recursive parameter estimation in a given interval was presented. Practical importance of this problem was stressed and simulation experiments support expected favourable influence on the transient identification behavior. Moreover at least the following interesting features of the problem assumed can be found: . It separates clearly the experience accumulation (computation of sufficient statistic) and decision making (selection of parameter estimates). . It shows that the information contained in the data must not be destroyed by any prior information. . It demands the use of the special type of mathematical programming in almost standard identification problem. This combination requires some modifications of the fundamental algorithm [ l ] , which result in still effective but more stable computation. As the by-product the re-computation of the LDLT decomposition after permuta-tions was solved in the unified and effective manner. This problem can be met also in computations connected with the system classification [9]. . It clarifies that the activity of any bound changes the other entries of point estimate. . It is useful also when computing time is bounded (this may appear in adaptive control), because of its monotonicity computation can be stopped before given number of iterations. At the end our main practical conclusion has to be repeated: any reasonable (may be heuristic) solution must respect the splitting of experience accumulation and decision making. (Received June 5, 1981.) R E F E R E N C E S [1] R. Fletcher, M. P. Jackson: Minimization of a quadratic function of many variables subject only to lower and upper bounds. J. Inst. Math. Appl. 14 (1976), 159—174. [2] V. Peterka: Bayesian approach to system identification. In: Trends and Progress in System Identification (P. Eykhoff, ed.), IFAC Series for Graduates, Research Workers and Practicing Engineers 1 (1981). [3] V. Peterka: Experience accumulation for decision making in multivariate time series. Problems Control Inform. Theory 7 (1978), 143-159. [4] M. H. De Groot: Optimal Statistical Decisions. McGraw-Hill, New York 1970. [5] M. Karny: Probabilistic Identification and Self-Reproducing Forms of Distribution Functions (in Czech). Ph. D. Dissertation, Institute of Information Theory and Automation of the Czechoslovak Academy of Sciences, Prague 1976. [6] L. Ljung: Analysis of vector stochastic algorithms. IEEE Trans. Autom. Control AC-22 (1977), 551-579. [7] G. K. Kel'mans, A. S. Poznyak, A. V. Chernitser: Adaptive locally optimal control. Internát. J. Systems Sci. 12 (1981), 235-254. [8] M. Kárný, J. Bohm, A. Halousková: Identification and Control of Cold Rolling System (in Czech). Internal Report o.p. Škoda Válcovny 1980. [9] M. Kárný: Automated structure determination of regression model (to appear). Ing. Miroslav Kárný ,CSc, Ústav teorie informace a automatizace ČSA V (Institute of Information Theory and Automation — Czechoslovak Academy of Sciences), Pod vodárenskou věži 4,182 08 Praha 8. Czechoslovakia. 49 

Recursive parameter estimation of regression model when the interval of possible values is given

http://dml.cz/bitstream/handle/10338.dmlcz/125590/Kybernetika_18-1982-1_4.pdf

Recursive parameter estimation of regression model when the interval of possible values is given

Abstract

Similar works

Full text

Available Versions

Institute of Mathematics AS CR, v. v. i.