Discovering Admissible Model Equations from Observed Data Based on Scale-Types and Identity Constrains

1 Introduction

The most well known pioneering system to discover sci­ entific law equations from experimental data is BA­ CON [Langley et al., 1985]. It searches for a complete equation governing the data measured in a continuous process, where the complete equation is an equation constraining n quantities with n - 1 degree of freedom 1 . FAHRENHEIT [Koehn and Zytkow, 1986], ABA­ CUS [Falkenhainer and Michalski, 1985], etc. are the successors that basically use similar algorithms to BA­ CON to discover a complete law equations. To re­ duce the high computational cost of their algorithm, some subsequent discovery systems, e.g., FAHREN­ HEIT, ABACUS and COPER [Kokar, 1985], introduced the use of the unit dimension of physical quantities to Yuji Niwa I.N.S.S., Inc. 64 Sata, Mihamacho, Mikatagun, Fukui 919-1205, Japan prune the meaningless solutions. A difficulty of this ap­ proach is the narrow applicability only to the quanti­ ties whose units are clearly known. On the other hand, the most recent scientific law discovery system, SDS, has overcome the difficulties of the past systems [Washio and Motoda, 1997] [Washio and Motoda, 1998]. It discov­ ers scientific law equations by limiting its search space to mathematically admissible equations in terms of the constraints of scale-type and identity. These constraints come from the basic characteristics of the quantities' def­ initions and the relations necessarily standing in the ob­ jective systems. The admissible equations discovered by SDS are considered to have valid structures reflecting the relations among quantities in the fundamental mecha­ nisms governing the objective system. The equations having such valid structures is called first principle equa­ tions in this paper. The detailed characterization of the first principle equations can be seen elsewhere [Washio and Motoda, 1998]. Since the knowledge of scale-types is widely obtained in various domains, SDS is applica­ ble to non-physical domains including biology, sociology, economics, etc.

1 The equation x 1 2 + x2 2 + ... + Xn 2 = 0 is not complete, since the values of all n quantities is 0, i.e., n quantities are constrained with no degree of freedom. On the other hand, X1 + X2 + ... + Xn = 0 is complete.

A major drawback of these approaches is the limited applicability to practical situations. They require the experimental environment and the interaction to control and measure the system states. The number of control­ lable quantities is quite limited, and even none of them are controllable due to some practical reasons in many scientific and engineering domains. For instance, the as­ tronomical experiments to control the parameters of fu­ sion reactions in the distant huge stars are physically impossible. The economical experiments to cause finan­ cial panics are unacceptable for our society. Under these situations where only passive observation is possible, the mathematical techniques, e.g., linear system identifica­ tion [Ljung, 1987] and neural network, have been tradi­ tionally applied to derive quantitative relations among observed quantities. However, the derived relations are not ensured to represent the first principle because they presume some structures of the model equations such as linear formulae and hierarchical sigmoid formulae. The discovery of the first principle equations under the pas­ sive observation will play highly important role to under­ stand the fundamental mechanisms underlying the variety of the objective systems. To achieve this aim by the technique of the aforementioned scientific discovery, the current framework must be changed to discover the first principle equations by using only the data obtained under the passive observation.

The past scientific discovery is for the class of the problem to discover the law equations under the exper­ imental environment. Its algorithm basically consists of two operations. The first is called bi-variate fitting which identifies the. relation within a pair of quantities, P;j = {x;,xj} ~X, where X= {xl,xz, ... ,xm} is the set of all quantities to represent the objective system. It derives the pairwise relation within P;j from the exper­ imental data in which the values of all quantities in the rest X - P;j is fixed by the experimental control. This pairwise relation is noted as fx-P;j(x;,xj) = 0. The bi­ variate fitting is required to identify the intrinsic struc­ ture of the relation within P;j under the exclusion of the influence from the other quantities. The second opera­ tion is to merge the multiple pairwise relations into an equation. Through the iteration of these two operations, the complete equation ¢>(x1 ,x2 , .·· ,xm) = 0 to represent the entire objective system is derived. In the new class of the problem to discover the law equations under the passive observation environment, the experimental con­ trol of the values of X- P;j is not allowed. Accordingly, the conventional bi-variate fitting is not applicable. In this paper, "quasi-bi-variate fitting", an extension of the bi-variate fitting based on a polynomial approximation, is proposed to enable the application of the framework of SDS to this class of the problem.

The proposed quasi-bi-variate fitting requires some as­ sumptions which are feasible in many practical applica­ tions. One is that the scale-types of all observed quanti­ ties are known. This does not limit the applicability of the proposed method because the scale-types of the mea­ surement quantities are widely known based on the mea­ surement theory as shown later. Another assumption is that the observed data are uniformly distributed over the value range that each quantity can take within the pos­ sible states of the objective system. If the observed data points are concentrated within the vicinity of a value for some quantity, the data set does not provide any meaningful information on the relation of the quantity with the others. Accordingly, the discovery of the first principle equations becomes difficult if this assumption is strongly violated. However, this requirements is not limited to our proposed approach. The lack of the uni­ form distribution of the data over a certain value range of a quantity implies the low observability of the quan­ tity [Ljung, 1987]. It is well known that the conventional approaches such as the linear system identification and the neural network do not derive valid models of the objective systems under the low observability condition. This limitation is generic for any data-driven modeling approaches, and further discussion on this issue is out of scope of this paper.

The objectives of this paper are (i) to propose the prin­ ciples and an algorithm of the quasi-bi-variate fitting under the framework of SDS, (ii) to evaluate the basic per­ formance of the proposed approach through simulations and (iii) to demonstrate its high practicality through a real application.

2 Background Principles

Before proposing quasi-bi-variate fitting, some back­ ground principles are explained to facilitate the compre­ hension. The details of the principles are described in our papers on SDS [Washio and Motoda, 1997] [Washio and Motoda, 1998]. Only its outline is explained in this section.

2.1 Scale-type Constraints

The rigorous definition of scale-type was given by Stevens [Stevens, 1946]. He defined the measurement process as "the assignment of numerals to object or events according to some rules." He claimed that dif­ ferent kinds of scale-types and different kinds of mea­ surement are derived if numerals can be assigned un­ der different rules, and categorized the scale-types of quantities based on the operation rule of the assign­ ment. Quantitative measurement quantities are mathe­ matically characterized and categorized into three major quantitative scale-types of interval scale, ratio scale and absolute scale. Examples of the interval scale quantities are temperature in Celsius and sound tone where the ori­ gins of their scales are not absolute, and are changeable by human's definitions. Its operation rule is "determina­ tion of equality of intervals or differences", and its ad­ missible unit conversion follows "Generic linear group: x' = kx + c". Examples of the ratio scale quantities are physical mass and absolute temperature where each has an absolute zero point. Its operation rule is "determina­ tion of equality of ratios", and its admissible unit con­ version follows "Similarity group: x' = kx". Examples of the absolute scale quantities are dimensionless quan­ tities. It follows the rule of "determination of equality of absolute value", and "Identity group: x' = x".

Luce claimed that the basic formula of the functional relation among quantities of ratio and interval scales can be determined by their scale-type features, if the quan­ tities are not coupled through any dimensionless quan­ tities [Luce, 1959]. Under this condition, the quantities should share some common basic dimensions, and con­ sequently the unit change of a quantity affects the value of other quantity. Suppose x and y are both ratio scale quantities, and y is defined by x through a continuous functional relation y = u(x). Suppose the form of u(x) is logarithmic, i.e., y = logx. We multiply a positive constant k to x, i.e., a change of unit, without violating the group structure of the ratio scale quantity x, then this leads u ( kx) = log k + log x. This fact causes the shift of the origin of y by log k, and violates the group structure of y which is the ratio scale quantity. Hence, the direct functional relation from x to y must not be logarithmic. Based on the admissibility condition of the relations among ratio and interval scale quantities, we mathematically derived the following two theorems to represent the generic formulae of the relations [Washio and Motoda, 1997]. 2

Theorem 1 (Extended Buckingham IT-theorem) If ¢(x1,x2,x3, ... ,xm) = 0 is a complete equation, and if each argument is one of interval, ratio and absolute scale-types, then the solution can be written in the form F(Ih,Ih, ... ,IIm-w) = 0,

where m is the number of arguments of¢, and w is the basic number of bases in x1, x 2, x3, ... , Xm, respectively.

Bases are such basic scaling factors and origins indepen­ dent of the other bases in the given ¢, for instance, as length [L], mass [M] and time [T] of physical unit and as temperature origin [to] of Celsius and elevation origin [ho] of potential energy for interval scale quantities. The relation of each IIi with the arguments of¢ is given by the following theorem.

Theorem 2 (Extended Product Theorem)

Assuming primary quantities in a set R are ratio scale­ type, and those in another set I are interval scale-type, the function p relating a secondary quantity II to Xi E R U I has one of the forms

where R and I can be null sets, P is a partition of I , and P 9 is a partition of I- I 9 where I 9 ~ I. All coefficients except II are constants.

The formula in Theorem 1 is called an "ensemble equa­ tion" and those in Theorem 2 "regime"s.

Table 1 shows all admissible bi-variate relations de­ duced from the "Extended Product Theorem". The co­ efficients Gij and Hij can be dependent on the other quantities except Xi and Xj. Thus, they are represented as Gij(X- Pij) and Hij(X- Pij), while aij is indepen­ dent, and remains constant. These consequences play an important role in the quasi-bi-variate fitting explained later.

2.2 Identity Constraint

When the scale-types of quantities are absolute and/or unknown as the case of "ensemble equation", the scale­ type constraints are not applicable. In such cases, the identity constraint is used to determine the admissible equation.

The basic principle of the identity constraints comes in by answering the question that "what kind of re­ lation holds among Bh , Bi and Bj , if Bi = ffli ( Bh) and Bj = f£J.(Bh) are known?'' For example, if Bi =

2 The original Buckingham IT-theorem [Buckingham, 1914] and Product Theorem [Bridgman, 1922] represent the generic relation among only ratio scale quantities.

LR is a set of pairwise terms having a bi-variate linear relation and LQ = U cELRC. P R is a set of pairwise terms having a bi-variate product relation and PQ = U cEPRC·

Ghi(Bj)Bh + Hhi(Bj) and Bj = Ghj(Bi)Bh + Hhj(Bi) are given, the following identity equation is obtained by solv­ ing each for Bh·

Because the third expression is linear with Bj for any Bi , the second must be so. Accordingly, the following must hold.

By substituting these to the second expression,

is obtained. Thus, by knowing some bi-variate linear relations among the quantities, the admissible equation formula for the whole quantities is derived.

This principle is generalized to various bi-variate rela­ tions f among multiple quantities. Table 2 shows such relations for linear relations and product relations.

3 Quasi-bi-variate Fitting

As noted in the first section, the conventional bi-variate fitting requires experimental control of some quantities, and is not applicable to the passive observation envi­ ronments. To overcome this difficulty, we propose the "quasi-bi-variate fitting" procedure which extracts a bi­ variate relation between two quantities under the ap­ proximated constant values of the other quantities.

3.1 Fitting for Scale-type Constraint

Figure 1 shows the outline of its principle for the ad­ missible bi-variate relations in Table 1. Let OBS = {X1 , X2 , ... , Xn} be a set of observations where each Xh(h = 1, ... , n) is a m-dimesional vector of observed values of the m quantities in X. The fitting of a can­ didate bi-variate formula for a pair of two quantities P;j = {x;,x·}(~ X) is applied to a subset of OBS. This subset DBS;jg is chosen in such a way that every quantity Xk E (X- P;j) takes a value in the vicinity of the value of Xkg, where X 9 = {xlg, Xzg, ... , Xmg} E OBS is an arbitrary chosen obser;vation vector. The vicinity of x kg is defined as

Ek determines the size of the vicinity. This vicinity is indicated by a rectangular cube in the upper figure of Fig. 1. Every admissible bi-variate formula indicated in Table 1 is generally represented by the form

Here, G;j and H;j are dependent on the quantities in X- P;j, while a;j remains constant. Given an OBS; · , if each Ek is moderately small, the values of G;j and R:j become slightly dependent on X- P;j, and their poly­ nomial approximation of the order p can be applied.

where G?jkg and Hi~kg stand for the coefficients of the h-th order of 6..xk at X 9 · The least square fitting of Eq.(3) approximately provides the functional relation within P;j and the coefficient a;j as depicted in the bot­ tom figure of Fig. 1. while almost excluding the influence of the other dimensions X - P;j.

After the least square fitting of this formula to OBSijg, the goodness of the fitting is evaluated by the following F -test. .

where

Here, u;/,u;; 2 and u00 2 are the correlation ofx; and Xj, the squared summation of x; and the squared summation of fitting error respectively. nijg is the total number of data points in OBSijg and F(1, nijg - 2, a) the lower bound ofF value under the degree of freedom (1, n;j 9 -2) and a risk rate a. The value of nijg is subject to the size of the vicinity EkS· a is set to be 0.05 throughout this paper. The quasi-bi-variate fitting and the F-test are repeated for multiple OBS;jJJS defined by q different X 9 s. This repetition is to confirm the stability of the F-test consequences. q is set to be 10 which is sufficient enough to check the stability of the consequences of the F-test. After these trials, the following x 2 -test over the q trials is conducted to check if a;j in Eq.(3) is identified as constant .

Here, 6aijg is the standard error of aijg estimated from the residual error of the quasi-bi-variate fitting. x 2 ( q -1, a) stands for the upper bound of x 2 value under the degree of freedom ( q -1) and the risk level a. The ap­ proximated bi-variate formula Eq.(3) that passed these tests is considered to be a part of the admissible model equation indicated in Theorem 1 and 2 . The expecta­ tion value of a;j is estimated as li;j = o:::::=l aijg)jp.

3.2 Fitting for Identity Constraint

For the bi-variate relations of the identity constraints indicated in the first column of Table 2, the similar scheme of the quasi-bi-variate fitting is applied. Let OBS = {01 , 0z, ... , 0n} be a set of observations, where each 0h(h = 1, ... , n) is a vector of observed values of them' quantities in 0 = {OI,Oz , ... ,Om'}· The fitting of a candidate bi-variate formula for a pair of two quanti­ ties P;i = {O;,Oi}(~ 0) is applied to OBS;ia· OBSijg is a subset of 0 B::; in the vicinity of the value of 0 kg, where 0 9 = {019 ,0z 9 , ... ,Om'g} is arbitrary chosen in OBS. Ev­ ery bi-variate formula indicated in Table 2 is generally represented by the form

Again, the terms of Gij(G-Pij) and Hij(8-Pij) are ap­ proximated by their polynomials. After the least square fitting of the approximated formula to OBSijg, the good­ ness of the fitting is evaluated by the F -test in the similar manner.

In the quasi-bi-variate fitting, the number of data in­ cluded in OBSijg increases by relaxing the size of the vicinity of X 9 and 8 9 · This has an effect of reducing the statistical error of the quasi-bi-variate fitting. On the other hand, if the size of the vicinity is too large, the higher order approximation is required to absorb the in­ fluence of the values of X- Pij and 8-P;j. However, ex­ cessively high order approximation may introduce some systematic error due to the over-fitting to the data. Ac­ cordingly, some appropriate values of Eks and p must be used for the given data.

4 Algorithm

As the details of the algorithm to discover a complete model equation in the frame work of SDS are represented in our previous paper, only its essential contents related to the quasi-bi-variate-fitting are explained in this sec­ tion [Washio and Motoda, 1997]. Initially, a set of ratio scale quantities RQ, a set of interval scale quantities IQ and a set of absolute scale quantities AQ which are re­ quired to express the objective model equation are given together with a set of observed data OBS of these quan­ tities.

Step (1-1)

The quasi-bi-variate fitting for scale-type constraints is applied to the bottom formula in Table 1 for pairwise interval scale quantities. The least-square fitting of the formula using the approximation of Eq.(3), F-test to check the goodness of fitting to the data of each subset OBSijg and x 2 -test to check the constant value of aij are conducted. Subsequently, the expectation value a;j is estimated, and the formulae together with the values of a;j are stored into an equation set I E.

This step is now demonstrated by an example of a moderately complex system depicted in Figure 2. This is an electric circuit where the model of this system based on the first principle is represented by the following equa­ tion involving eight quantities.

where RBE is the resistance between the base and the emitter and hp E the ratio between the base current and the collector current, respectively. Vi and Vz are in­ terval scale, and hp E is absolute scale. The rests are ratio scale. The observed data set is obtained by a numerical simulator. The values of parameter quan­ tities are set to be R2 = 1000S1, RBE = 10 6 S1 and hFE = 100. The value ranges of the variable quantities are taken to be OQ < RI < lOOOQ, OQ < R3 < 1000Q

and OV < V2 < VI < 30V. The values of these vari­ ables are generated by using uniform random numbers over their value ranges in the simulation. Only the five variable quantities RI, R3 , VI, V2 and I are assumed ob­ servable in this demonstration. Thus, IQ = {VI, V2}, RQ = { RI, R3 , I} and AQ = ¢. The values of the pa­ rameter quantities are implicitly assumed to be constant. The total number of data points provided in OBS is 500, and no observation noise is added here. Our proposed method has been implemented in a prototype program. The size of each vicinity Ek, has been set at 15% of the difference between the maximum and the minimum val­ ues of Xk. The Oth order approximation is used because the observation is not distorted by any noise in this case. As the interval scale quantities are limited to Vi and V2, the bottom linear equation in Table 1 is immediately ap­ plied to this pair, and the relation 0.98VI - Vz = Gv, v 2 is identified. Thus IE= {0.98VI - V2 = Gv, v 2 }.

Step (1-2)

This step first applies the following triplet-tests. For a triplet of the linear formulae among {xi, Xj , xh} in IE,

if they are mutually consistent in terms of as, the fol­ lowing condition should be met.

Because of the existence of the noise and the fitting error, this condition does not hold in exact manner, even if the three formulae are consistent. Thus, the following normal distribution test judges if the l.h.s. and the r.h.s. of the above expr~ssion are equal.

where

Here, 8aij = JC2:,~= 1 6a~i 9 )/p, and 8ajh and 8ahj are similarly defined. N(O, a 2 , a/2) stands for the upper bound of the error under the normal distribution and the risk level a. This test is applied to every triplet of equations in IE, and every maximal convex set MCS is searched. A convex set is a set where each triplet of equations among the quantities in this set has passed the test Eq.(8). And, the maximal convex set MCS is a con­ vex set where any superset of the set is not a convex set. In addition, every formula in IE which does not belong to any consistent triplet is also regarded as a tiny MCS. Once all MCSs are found, the formulae are merged into the following form in every M C S.

where r is an intermediate quantity which appears in the reasoning process. Before the final value of a, is de­ termined, the following integer-test is applied.

This is based on the observation that the majority of the first principle based equations have integer power co­ efficients and integer linear coefficients for interval scale quantities.

In the current circuit example, an MCS is uniquely determined because IE= {0.98V1 - V2 = Gv,v 2 } con­ tains only one formula, and thus the triplet test is not required. The example of the triplet test is shown in section 6. The above integer-test set av,, the coefficient of V1 , to be 1 because 28av, = 0.092. Thus, we obtain IE= {V1 - V2 = Gv, v 2 }. Furthermore, V1 and V2 in IQ is merged into Gv, v,, and Gv, v 2 is stored into a quan­ tity set TQ as TQ = { Gv, v 2 }. Gv, v, is a new ratio scale quantity by the mutual cancellation of the basic origins of V1 and V2. Finally, TQ = TQ + RQ becomes {R1, R3,I , Gv, v 2 }.

Step (2-1)

Similarly to step (1- 1), the quasi-bi-variate fitting, F­ test and x 2 -test are performed on the quantities included in TQ. Then, the discovered equations are stored in an equation set RE. The unique difference from step (1- 1) is to apply the formulae except the bottom one in Table 1 in the quasi-bi-variate fitting. In the circuit example, only the pair of I and Gv, v 2 is found to satisfy the first formula in Table 1 as I = G rev, v 2 �k 0 ~;. Thus, RE = {I= G IGv, v2 �~,o~n.

Step (2-2)

The triplet-test among the formulae in RE is con­ ducted. The basic procedure is identical with the step (1-2). In the example, as RE contains only one for­ mula again, a unique MCS becomes {I,Gv 1 v2!, and they are merged into a term G rev, v 2 = I /�.;, 0 ~ 7 · The value acv, v 2 = 1.003 is modified to 1 in the mteger test since 28acv 1 v 2 = 0.014. Thus, G IGv 1 v 2 = I/ Gv, v 2 · Then, TQ becomes {R1, R3 , G rev, v 2 }, and finally TQ = TQ+AQ = {Rl , R3,Grcv, v 2 }·

Step (3)

This is the step to apply the quasi-bi-variate fitting for identity constraints. A formula is arbitrary selected from the first column of Table 2. In our current program, a linear formula Bj = G;jBi + H;j has the first priority in the selection and a power formula Bj = H;jB~'j the next priority. The quasi-bi-variate fitting of the for­ mula and F-test are applied to each pair of quantities in TQ. If some pairs of the quantities are judged to well fit to the bi-variate formula, the identity constraints is applied. In the example of TQ = { R1 , R3, G rev, vJ, the bi-variate linear relations in the pairs of { R1 , G IGv v } 1 2 and { R3, G IGv 1 v 2 } are accepted through the F -test, and thus the following multi-linear formula obtained from the principle of the identity constraints is applied to the en­ tire data set.

Its least square fitting is accepted by F-test, and thus AE = {0 = ao + a1 G IGv 1 v 2 + a2R1 + a3R3 + R1R3}. The values ao = 10 7 , a1 = -10 3 , a2 = 10 4 and a3 = 10 3 are obtained by the integer-test. By substituting the formulae in IE and RE to G rev, v 2 , the final solution of the admissible model equation of

is resulted. As th values of a0 - a3 correspond to the relations ao = RBER2/hFE, a1 = -R2, a2 = RBE/hFE and a3 = R2 , this equation is known to be equivalent to Eq.(7).

Step (4)

Finally, when multiple candidate model equations re­ main, a parsimony criterion is applied to prioritize the candidates. Though MDL principle is a representative criterion, AIC, which is widely used in statistics to de­ termine an appropriate numerical model equation, is ap­ plied in our program [Akaike, 1978]. The index of AIC is calculated through the expression

where n = IOBSI, Ve the residual error variance of the model equation and M the number of the coefficients included in the model. The model equation having less value of AIC is preferred in the sense of the parsimony criterion. This is not used in the example of the cir­ cuit, since the unique solution Eq.(10) is obtained. Its example is given in latter· section 6.

5 Evaluation through Simulation

Table 3 indicates the required computation time for var­ ious examples. "Ideal Gas" is the simulation of the state equation of the ideal gas. "Coulomb", "Stoke's" and "Momentum" are the simulations of Coulomb force law, Stoke's equation and the momentum balance equation. "Circuit*1" is the case of the aforementioned electric cir­ cuit where R2 , RBE and hFE are hidden parameters, and "Circuit*2" is the case of the identical circuit where all quantities are observable. They are represented by vari­ ous number of quantities. The computation time of the proposed algorithm has been evaluated for various num­ bers of the data points for each example. The computa­ tion time does not change significantly with the increase of the data size. This is because the 10 vicinities selected in the quasi-bi-variate fitting cover only a limited por­ tion of the given data when the data size is large. Thus the required computation time increases very slowly. In contrast, the computation time is sensitive to the size of the objective system. The increase is almost order of O(m 2 ) where m is the number of the quantities in the data. This is due to the quasi-bi-variate fitting hav­ ing the complexity O(m 2 ). Though the most complex process is the triplet-test which is O(m 3 ), this test is very simple compared with the data fitting to many data points.

Table 4 shows the average of the relative error of the coefficients under several conditions of the data size, the relative noise level and the order of the approximation of the quasi-bi-variate fitting in case of the aforemen­ tioned electric circuit. The size of the vicinity Ek is kept at 15%. When the amount of the data is very limited, the error rate increases significantly. This is because the data points covered by a vicinity is so small that the suffi­ cient statistic accuracy is not maintained. The accuracy of the coefficients is also influenced by the approximation order. In general, the 1st order approximation shows the good performance. This tendency becomes significant, when the number of the data is very limited, and/or the noise level is high. The Oth order approximation does not effectively reduce the influence of (X - Pij ), when the data points are sparse in the state space. Besides, the 2nd order approximation also becomes erroneous due to the over fitting effect, when the data is sparse and/ or the noise level is high.

6 Application to Practical Problem

The proposed method has been applied to a real world problem. The objective of the application is to discover a generic law formula governing the mental preference of people on their houses subject to the cost to buy the house and the social risk at the place of the house. We designed a questionnaire sheet to ask the preference of the house in the trade off between the frequency of huge earthquakes, x 1 (earthquake/year), and the cost to buy, x2 ($). In the questionnaire, 9 combinations of the cost and the earthquake frequency are presented, and each person chooses its preference from the 7 levels for each combination. We distributed this questionnaire sheet to the people owning their houses in the s1,1burb area of Tokyo, and totally 400 answer sheets are collected back. The answer data has been processed by following the method of successive categories which is widely used in the experimental psychology to compose an interval scale preference index y [Torgerson, 1958], and OBS = {X1,X2, ... , X4oo} where Xi= [xli,X2i,Yi] is obtained.

The proposed method has been applied to figure out a law equation y = f(x 1 ,x2 ), where x1 and x2 are ra­ tio scale quantities. Here, RQ = {x1 ,x2 },IQ = {y} and AQ = </J . Because IQ contains only one quantity, steps (1- 1) and (1- 2) are skipped, and the quasi-bi­ variate fitting of 1st order approximation is applied to RQ = RQ + IQ = {x1 , x2 , y} in step (2- 1). First, the top formula in Table 1 is tested for the relation between x1 and x2 , and x1 = a(y )x; 0 · 25 is obtained. Next, the second and third formulae are tested for x1 andy. Then, y = a(x2 )x- 0 · 23 + b(x2 ) andy= 0.62log x 1 + b(x2 ) have been identiRed. Similar search has been made for x2 and y, andy = a(x1 )x~· 026 +b(xl) andy = 0.34log x2 + b(x1 ) are derived. In step (2 - 2), the triplet-test among { x1 , x2 , y} is conducted. As Theorem 2 limits the ad­ missible formulae for these quantities to y = bxr' x~ 2 + c and y = a1 log x1 + a2 log x2 + c, the consistency among the coefficients obtained in step (2 - 1) are checked for these formulae. As the result , the consistency has been confirmed for both. Consequently, we obtained the fol­ lowing two candidates.

Steps (3- 1) and (3- 2) are skipped, since the other quantities to be merged do not exist. In step ( 4), the values of AIC are evaluated for each candidate, where the former Eq.(12) is preferred. Eq.(12) can evaluate the subjective preference in the accuracy of almost ±1 preference levels of the questionnaire from the values of x1 and x2.

7 Discussion and Related Work

A scientific discovery system called LAGRANGE [Dze­ roski and Todorovski, 1994] is also applicable to the con­ dition of the passive observation. It uses the principles of ILP and generate/test. Though no equation classes are presumed in this approach, many spurious solutions can be derived due to the weakness of the search heuris­ tics. Also, it indicates high computational complexity. TETRAD [Glymour, 1995] is another system to iden­ tify the models of the objective system from the passive observation. Its basic framework takes the bottom up modeling approach. However, the class of the model formulae is presumed such as linear expressions. In con­ trast, the method proposed in this paper has a strong mathematical background to characterize first principle equations. Moreover, it has a high applicability to the passive observation data, while maintaining the flexibil­ ity of the bottom up modeling approach taken by the conventional scientific discovery systems.

The source of the advantage of our proposed method is the systematic use of the constraints of scale-types and identity with the approximation in quasi-bi-variate fitting. This is considered to be a typical example that Ginsberg claimed [Ginsberg and Geddis, 1991]. He claimed that any domain-dependent control rules can be replaced with a domain-independent control rules and modal sentences describing the structure of the search space. The knowledge of the scale-types and the quasi­ bi-variate approximation have been implicitly used by scientists as domain-dependent control rules of their rea­ soning. In our work, these rules have been replaced as Ginsberg claimed. The constraints and the approxima­ tion are formalized as generic domain-independent con­ trol rules applicable to any objective system represented by numerical quantities. The modal knowledge required to control the reasoning by these generic rules is con­ centrated on the scale-type information of each quan­ tity and the empirical quasi-bi-variate relation. On the other hand, Minton argued that in many cases, domain­ dependent control rules cannot, in a practical sense, be derived due to the complexity of the reasoning that would be required [Minton, 1996]. Since the concepts such as the scale-types and the quasi-bi-variate approx­ imation have been established based on massive experi­ ence of the scientists for hundreds of years, his argument also holds.

8 Conclusion

In this paper, quasi-bi-variate fitting, an extension of the bi-variate fitting based on an polynomial approximation, has been proposed. This extension enables to handle the new class of the problem to discover the law equations under the passive observation environment. Its basic performances in terms of computation time and noise ro­ bustness have been evaluated through simulations. The evaluation indicates the satisfactory performance to dis­ cover the model equation based on the first principle of objective system of moderately large size under practical noise levels. Finally, a real application to discover a law equation in socio-psychology was demonstrated, and its practicality has been readily confirmed.