The present invention claims priority to Japanese Patent Application No. 2014-207163 filed Oct. 8, 2014, which is incorporated herein by reference in its entirety.
1. Technical Field
The present invention relates to a technique of generating a calibration curve used for deriving the content of a target component related to a subject from observation data of the subject and a technique of acquiring the content of the target component related to the subject.
2. Related Art
In the related art, the use of independent component analysis for generating a calibration curve of a target component and calibration of a target component using the calibration curve has been suggested. For example, in Japanese Patent No. JP-A-2013-160574, the independent component analysis related to observation data of a plurality of samples whose components are known is performed at the time of generating a calibration curve and a simple regression calibration curve is determined from a relationship between the mixing ratio of acquired independent components and the amount of known target components. In addition, at the time of calibration of the target component, calibration of target components is performed using the independent components and the calibration curve acquired at the time of generating a calibration curve with respect to the observation data related to a target sample whose component amount is unknown. In addition, a technique of removing a fluctuation in the observation data by applying project on null space (PNS) or standard normal variate transformation (SNV) and performing calibration with high precision as first pre-processing of the observation data is described.
Despite these advances, however, there is problem with the known methods and systems because fluctuation cannot be completely removed and the precision cannot be sufficiently obtained in some cases in the pre-processing of the related art.
An advantage of some aspects of the invention is to solve at least a part of the problems described above, and the invention can be implemented as the following aspects or application examples.
(1) A first aspect of the invention provides a calibration curve generation method of generating a calibration curve used for deriving a content of a target component of a subject from observation data of the subject. The calibration curve generation method includes acquiring an independent component matrix including independent components of each sample, by performing first pre-processing that includes normalizing the observation data, second pre-processing that includes whitening, and independent component analysis processing in order. Further, the same noise is added to the observation data related to the plurality of samples in the first pre-processing.
According to this method, by adding the same noise to the observation data related to the plurality of samples, it is possible to reduce influence of the fluctuation in the observation data compared to a case where noise is not added thereto and to improve the calibration precision.
In one embodiment, the calibration curve generating method includes (a) acquiring, by a computer, the observation data related to the plurality of samples of the subject; (b) acquiring, by the computer, the content of the target component related to each of the samples; (c) estimating, by the computer, a plurality of independent components when the observation data for each of the samples is divided into the plurality of independent components and acquiring a mixing coefficient corresponding to the target component for each of the samples based on the plurality of independent components; and (d) acquiring, by the computer, a regression formula of the calibration curve based on the content of the target component of the plurality of samples and the mixing coefficient for each of the samples. The (c) estimating of the plurality of independent components includes (i) acquiring, by the computer, an independent component matrix that includes the independent components of each of the samples; (ii) acquiring, by the computer, an estimated mixed matrix indicating a set of vectors that define a ratio of independent component elements for each of the independent components in each of the samples from the independent component matrix; (iii) acquiring, by the computer, a correlation with respect to the content of the target component of the plurality of samples for each of the vectors included in the estimated mixed matrix and selecting the vector determined that the correlation is the maximum as a mixing coefficient corresponding to the target component. The independent component matrix is acquired by the computer performing first pre-processing that includes normalizing the observation data, second pre-processing that includes whitening, and independent component analysis processing in order, and the computer adds the same noise to the observation data related to the plurality of samples in the first pre-processing.
According to this calibration curve generation method, a calibration curve for deriving the amount of target components included in the subject from the observation data of the subject is generated from the observation data acquired from each of the samples and the content of the target component related to the plurality of samples of the subject. For this reason, when the calibration curve is used, the content of the target component can be acquired with high precision even if the observation data of the subject is one piece of data. Accordingly, when the calibration curve is generated in advance according to the calibration curve generation method, only one piece of observation data related to the subject needs to be acquired during the calibration. As a result, the amount of target components can be acquired with high precision from one piece of observation data which is a measured value. In addition, since an estimated mixed matrix is acquired and then a vector highly correlative with the content of the target components of samples is drawn out from the estimated mixed matrix, a mixing coefficient with high estimation precision can be obtained. Moreover, the calibration precision can be improved by adding the same noise to the observation data related to the plurality of samples in the first pre-processing. Further, the computer described in the processes (a) to (d) may be the one same computer or a plurality of computers different from one another. That is, each of the processes may be performed through data communication of the plurality of computers. Similarly, when each of the processes is performed by a computer in the present specification, the computer may be the one same computer or a plurality of computers different from one another.
(2) A second aspect of the invention provides a calibration curve generation device which generates a calibration curve used for deriving a content of a target component of a subject from observation data of the subject. The calibration curve generation device includes an independent component matrix calculation unit that acquires an independent component matrix having independent components of each sample. The independent component matrix calculation unit acquires the independent component matrix by performing first pre-processing that includes normalizing the observation data, second pre-processing that includes whitening, and independent component analysis processing in order. Further, the independent component matrix calculation unit adds the same noise to the observation data related to the plurality of samples in the first pre-processing.
According to this calibration curve generation device, by adding the same noise to the observation data related to the plurality of samples, it is possible to reduce influence of the fluctuation in the observation data compared to a case where noise is not added thereto and to improve the calibration precision.
In one embodiment, the calibration curve generation device includes a sample observation data acquisition unit that acquires the observation data related to a plurality of samples of the subject, a sample target component amount acquisition unit that acquires the content of the target component related to each of the samples, a mixing coefficient estimation unit that estimates a plurality of independent components when the observation data for each of the samples is divided into the plurality of independent components and acquires a mixing coefficient corresponding to the target component for each of the samples based on the plurality of independent components, and a regression formula calculation unit that acquires a regression formula of the calibration curve based on the content of the target component of the plurality of samples and the mixing coefficient for each of the samples. The mixing coefficient estimation unit includes an independent component matrix calculation unit that acquires an independent component matrix including each of the independent components of each of the samples, an estimated mixed matrix calculation unit that acquires an estimated mixed matrix indicating a set of vectors that define a ratio of independent component elements for each of the independent components in each of the samples from the independent component matrix, and a mixing coefficient selection unit that acquires a correlation with respect to the content of the target component of the plurality of samples for each of the vectors included in the estimated mixed matrix and selecting the vector determined that the correlation is the maximum as a mixing coefficient corresponding to the target component. The independent component matrix calculation unit acquires the independent component matrix by performing first pre-processing that includes normalizing the observation data, second pre-processing that includes whitening, and independent component analysis processing in order. Further, the independent component matrix calculation unit adds the same noise to the observation data related to the plurality of samples in the first pre-processing.
According to this calibration curve generation device, a calibration curve for deriving the amount of target components included in the subject from the observation data of the subject is generated from the observation data acquired from each of the samples and the content of the target component related to the plurality of samples of the subject. For this reason, when the calibration curve is used, the content of the target component can be acquired with high precision even if the observation data of the subject is one piece of data. Accordingly, when the calibration curve is generated in advance according to the calibration curve generation method, only one piece of observation data related to the subject needs to be acquired during the calibration. As a result, the amount of target components can be acquired with high precision from one piece of observation data which is a measured value. In addition, since an estimated mixed matrix is acquired and then a vector highly correlative with the content of the target components of samples is drawn out from the estimated mixed matrix, a mixing coefficient with high estimation precision can be obtained. Moreover, the calibration precision can be improved by adding the same noise to the observation data related to the plurality of samples in the first pre-processing.
(3) A third aspect of the invention provides a target component calibration method of acquiring a content of a target component related to a subject. The target component calibration method includes acquiring a mixing coefficient with respect to the target component related to the subject based on observation data and calibration data related to the subject, first pre-processing that includes normalizing the observation data and second pre-processing that includes whitening are performed in order in the acquiring of the mixing coefficient, and noise is added to the observation data related to the subject in the first pre-processing.
According to this target component calibration method, by adding the noise to the observation data related to the subject, it is possible to reduce influence of the fluctuation in the observation data compared to a case where noise is not added thereto and to improve the calibration precision.
In one embodiment, a target component calibration method includes (a) acquiring, by a computer, observation data related to the subject; (b) acquiring, by the computer, calibration data including at least an independent component corresponding to the target component; (c) acquiring, by the computer, a mixing coefficient corresponding to the target component related to the subject based on the observation data and the calibration data related to the subject; and (d) calculating, by the computer, the content of the target component based on a constant of a regression formula indicating a relationship between the mixing coefficient and the content corresponding to the target component, prepared in advance, and the mixing coefficient acquired in (c) the acquiring of the mixing coefficient. The computer performs first pre-processing that includes normalizing the observation data and second pre-processing that includes whitening in this order in (c) the acquiring of the mixing coefficient, and the computer adds noise to the observation data related to the subject in the first pre-processing.
According to this target component calibration method, the content of the target component related to the subject can be acquired with high precision only by obtaining one piece of observation data related to the subject. Further, by adding noise to the observation data related to the subject, the calibration precision can be improved.
(4) A fourth aspect of the invention provides a target component calibration device which acquires a content of a target component related to a subject. The target component calibration device includes a mixing coefficient calculation unit that acquires a mixing coefficient with respect to the target component related to the subject. The mixing coefficient calculation unit performs first pre-processing that includes normalizing observation data and second pre-processing that includes whitening in order and adds noise to the observation data related to the subject in the first pre-processing.
According to this target component calibration device, by adding the noise to the observation data related to the subject, it is possible to reduce influence of the fluctuation in the observation data compared to a case where noise is not added thereto and to improve the calibration precision.
In one embodiment, a target component calibration device includes a subject observation data acquisition unit that acquires the observation data related to the subject, a calibration data acquisition unit that acquires calibration data including at least an independent component corresponding to the target component, a mixing coefficient calculation unit that acquires a mixing coefficient corresponding to the target component related to the subject based on the observation data and the calibration data related to the subject, and a target component amount calculation unit that calculates the content of the target component based on a constant of a regression formula indicating a relationship between the mixing coefficient and the content corresponding to the target component, prepared in advance, and the mixing coefficient acquired by the mixing coefficient calculation unit. The mixing coefficient calculation unit performs first pre-processing that includes normalizing the observation data and second pre-processing that includes whitening in order and adds noise to the observation data related to the subject in the first pre-processing.
According to this target component calibration device, the content of the target component related to the subject can be acquired with high precision only by obtaining one piece of observation data related to the subject. Further, by adding noise to the observation data related to the subject, the calibration precision can be improved.
In the method or device, the process using the project on null space and the normalization may be performed in order after the noise is added in the first pre-processing. According to this configuration, the influence of various fluctuations in the observation data can be effectively reduced and the calibration precision can be improved by both effects of the process of adding the noise and the process using the project on null space.
Alternatively, in the method or device, the addition of the noise and the normalization may be performed in order after the process using the project on null space is performed in the first pre-processing. According to this configuration, the influence of various fluctuations in the observation data can be effectively reduced and the calibration precision can be improved by both effects of the process of adding the noise and the process using the project on null space.
The invention can be implemented using various aspects other than the above-described aspects. For example, the invention can be implemented using an electronic device including the above-described device, a computer program implementing functions of respective units of the above-described device, and a recording medium which stores a computer program and is not temporary (non-transitory storage medium).
The invention will be described with reference to the accompanying drawings, wherein like numbers reference like elements.
Hereinafter, embodiments of the invention will be described in the following order.
A. Outline of calibration curve generation process and calibration process:
B. Calibration method using noise addition:
C. Calibration curve generation method:
D. Calibration method of target component:
E. Various algorithms and influence thereof:
F. Modification Examples:
In the present embodiment, the following abbreviations are used.
When a calibration curve is generated, first, fluctuation or noise included in observation data is reduced by performing pre-processing on the observation data (
Next, as shown in (D) to (F) in
The simple regression formula of a calibration curve is represented by the following expression.
C=u·P+v (1)
Here, C represents the content of a target component, P represents the inner product (mixing coefficient), and u and v represent a constant.
In the present embodiment, as described below, the influence due to various fluctuations in measurement data (observation data) is reduced and the calibration precision is improved by adding noise to the measurement data in at least one of the pre-processing of
In the calibration curve generation process shown in
For example, the noise addition is performed as follows. A noise vector Z according to normal distribution with an average μ (=0), a distribution σ2 (predetermined set value), and a data length N is generated using a computer program having a function of generating random numbers or a function. Here, N represents the data length (number of segments in the wavelength bandwidth) of measurement data. The term “noise vector Z” indicates one-dimensional data having N number of elements. Further, the value of the distribution σ2 of the noise vector Z may be experimentally adjusted along with the kind of sample or the waveform of measurement data. It is expected that the calibration precision is further improved through the adjustment.
Next, a fixed noise matrix Nz is generated using the noise vector Z.
N
z
=Z×1T (2)
Here, 1T represents a column vector in which n elements are all 1.0 and n represents the number of pieces of measurement data (that is, the number of samples). Accordingly, the fixed noise matrix Nz is a matrix with n rows and N columns, in which row vectors are all the same noise vector Z.
Next, a measurement data matrix X′ with noise is generated by adding the fixed noise matrix Nz to a measurement data matrix X.
X′=X+N
z (3)
Here, the measurement data matrix X is a matrix with n rows and N columns in which each of n pieces of measurement data is set as a row vector. In the calibration curve generation procedures described below, the measurement data matrix X or the measurement data matrix X′ with noise is referred to as “spectral data X.” Another pre-processing or independent component analysis described in (A) to (C) in
In the calibration process shown in (A) to (D) in
In the first sample group SG1, the calibration precision SEP is 32.75 mg/dL in a case where noise is not added, but improved to 21.13 mg/dL in a case where noise is added. Further, the correlation coefficient R2 is improved to 0.9814 from 0.9576. In the second sample group SG2, the calibration precision SEP is 22.01 mg/dL in a case where noise is not added, but improved to 14.14 mg/dL in a case where noise is added. Further, the correlation coefficient R2 is improved to 0.9916 from 0.9799. In the third sample group SG3, the calibration precision SEP is 22.13 mg/dL in a case where noise is not added, but improved to 11.33 mg/dL in a case where noise is added. Further, the correlation coefficient R2 is improved to 0.9946 from 0.9798.
In this manner, it can be confirmed that the calibration precision is improved by adding noise to the measurement data in the pre-processing.
The process 1 is a preparation process and performed by an operator. The operator prepares the same kind of a plurality of samples (for example, a glucose aqueous solution or a human body) in which the contents of target components are different from one other. In the example, n (n represents an integer of 2 or greater) samples are used.
Process 2
The process 2 is a process of measuring the spectrum and performed by the operator using a spectrometer. The operator measures the spectrum of spectral reflectance related to each of the samples by imaging each of the plurality of samples prepared in the process 1 using the spectrometer. The spectrometer is a known device in which light from an object to be measured passes through a spectroscope, the spectrum output from the spectroscope is received by an imaging surface of an imaging element, and thus the spectrum is measured. The spectrum of the spectral reflectance and the spectrum of the absorbance satisfy the relationship represented by the following expression.
[Absorbance]=−log10[Reflectance] (4)
The spectrum of the measured spectral reflectance is transformed to the absorbance spectrum using the expression (4). The reason for conversion of the spectrum of the spectral reflectance to the absorbance spectrum is that a linear combination needs to be established in a mixed signal to be analyzed by the independent component analysis described below. In this instance, the linear combination related to the absorbance is established according to a Lambert-Beer's law. Therefore, in the process 2, the absorbance spectrum may be measured in place of the spectral reflectance spectrum. As the measurement results, absorbance distribution data showing characteristics with respect to the wavelength of the object to be measured is output. The absorbance distribution data is referred to as spectrum data.
Moreover, the spectral reflectance spectrum or the absorbance spectrum may be estimated from other measurement values instead of measuring these spectra using a spectroscope. For example, a sample is measured using a multi-band camera and the spectral reflectance or the absorbance spectrum may be estimated from an obtained multi-band image. As such an estimation method, a method described in JP-A-2001-99710 can be used.
The process 3 is a process of measuring the content of target components and is performed by the operator. The operator performs chemical analysis on each of the plurality of samples prepared in the process 1 and measures the content of the target components (for example, the amount of glucose) in regard to each of the samples. In a case where the content of the target components in the samples prepared in the process 1 is known, the process 3 can be omitted.
The process 4 is a process of estimating a mixing coefficient and is typically performed using a computer.
The computer 100 is a known device including a CPU 10 that performs various processes or control by executing computer programs, a memory 20 (storage unit) which stores a location of data, a hard disk drive 30 that stores computer programs and data, an input I/F 50, and an output I/F 60.
The spectrometer 200 shown in
As a result of acquisition of the spectrum data and the content of target components described above, a dataset (hereinafter, referred to as a “measurement dataset”) DS1 including the spectrum data and the content of target components is stored in the hard disk drive 30 of the computer 100.
The CPU 10 performs a process of estimating a mixing coefficient, which is the operation of the process 4, by loading a predetermined program stored in the hard disk drive 30 in the memory 20 and executing the program. Here, the predetermined program can be downloaded using a network such as the Internet. In the process 4, the CPU 10 functions as the mixing coefficient estimation unit 430 of
The independent component analysis (ICA) is one of multi-dimensional signal analysis methods and is a technique of observing a mixed signal in which independent signals overlap each other under several different conditions and separating independent original signals based on the observation. When the independent component analysis is used, the spectrum of the independent component can be estimated from the spectrum data (observation data) obtained in the process 2 by grasping the spectrum data obtained in the process 2 as data mixed with m independent components (unknown) including target components.
In the present embodiment, the independent component analysis is performed by the three processing units 450, 460, and 470 shown in
In addition, in a case where the SNV 452 is performed on the spectrum data obtained in the process 2 of
Moreover, as the first pre-processing, a process other than the noise addition, the SNV, and the PNS may be performed. Preferably, any normalization process is performed in the first pre-processing, but the normalization process may be omitted. Hereinafter, the first pre-processing unit 450 is referred to as the “normalization processing unit.” The contents of these two processes 452 and 454 will be described below. Further, in a case where the data to be processed provided for the independent component matrix calculation unit 432 is normalized data, the first pre-processing may be omitted.
In the second pre-processing unit 460, pre-processing using any one of principal component analysis (PCA) 462 and factor analysis (FA) 464 can be performed. In addition, as the second pre-processing, processes other than PCA or FA may be performed. Hereinafter, the second pre-processing unit 460 is referred to as the “whitening processing unit.” In a general technique of the ICA, dimension compression of data to be processed or decorrelation is performed as the second pre-processing. Since a transformation matrix to be acquired by the ICA is restricted to an orthogonal transformation matrix by the second pre-processing, the calculation amount of the ICA can be reduced. Such second pre-processing is referred to as “whitening,” and the PCA is used in many cases. However, in a case where random noise is included in the data to be processed, an error may be generated in the results of the PCA due to the influence of the random noise. Here, in order to reduce the influence of the random noise, it is preferable to perform whitening using the FA having robustness with respect to noise in place of the PCA. The second pre-processing unit 460 of
The independent component analysis processing unit (ICA processing unit) 470 estimates the spectrum of independent components by performing the ICA with respect to the spectrum data to which the first pre-processing and the second pre-processing are applied. The ICA processing unit 470 can perform analysis using any one of a first processing 472 using a kurtosis as an independence index and a second process 474 using β divergence as the independence index. As the index for separating independent components from each other, the ICA uses higher order statistics representing independence of separated data as the independence index. The kurtosis is a typical independence index. However, in a case where an outlier such as spike noise is present in the data to be processed, the statistics including the outlier are calculated as the independence index. For this reason, an error is generated between original statistics related to the data to be processed and calculated statistics and this causes degradation of separation precision. Here, in order to reduce the influence from the outlier in the data to be processed, it is preferable to use an independence index which is unlikely to be affected by the outlier. It is possible to use the β divergence as the independence index having such characteristics. The contents of the kurtosis and the β divergence will be described below. Further, as the independence index of the ICA, an index other than the kurtosis and the β divergence may be used.
The contents of a typical process of the independent component analysis will be described below. It is assumed that a spectrum S (hereinafter, also simply referred to as an “unknown component”) of m unknown components (sources) is provided as a vector of the expression (5) below and n pieces of spectrum data X obtained in the process 2 is provided as a vector in the expression (6) below. Further, respective elements (S1, S2, . . . , Sm) included in the expression (5) are vectors (spectra). That is, the element S1 is represented by the expression (7). Elements (X1, X2, . . . , Xn) included in the expression (6) are vectors and, for example, the element Xj is represented by the expression (8). The suffix j of the element Xj represents the number of wavelength bandwidths measuring a spectrum. Moreover, an element number m of the spectrum S of an unknown component represents an integer of 1 or greater and is empirically or experimentally determined in advance according to the kind of sample.
S=[S1, S2, . . . , Sn]T (5)
X=[X1, S2, . . . , Xn]T (6)
S1={S11, S12, . . . , S11} (7)
X1={X11, X12, . . . , S11} (8)
The respective unknown components are statistically independent. The unknown component S and the spectrum data X satisfy a relationship of the following expression.
X=A·S (9)
In the expression (9), A represents a mixed matrix and can be represented by the expression (10) below. Further, the character “A” here needs to be written in bold as shown in the expression (10), but a normal character is used here because of the restriction of characters used in the specification. Hereinafter, similarly, other bold characters indicating matrixes are written in normal characters.
A mixing coefficient aij included in the mixed matrix A indicates a degree of contribution of an unknown component Sj (j=1 to m) to spectrum data Xi (i=1 to n) which is observation data.
In a case where the mixed matrix A is known, a least square solution of the unknown component S can be simply acquired as A+·X using a pseudo inverse matrix A+ of A. However, in the present embodiment, since the mixed matrix A is unknown, the unknown component S and the mixed matrix A need to be estimated from only the observation data X. That is, as represented by the expression (11) below, a matrix (hereinafter, referred to as an “independent component matrix”) Y indicating a spectrum of independent components is calculated using a separation matrix W of m×n from only the observation data X. As an algorithm acquiring the separation matrix W in the expression (11) below, various algorithms such as Infomax, fast independent component analysis (FastICA), and joint approximate diagonalization of eigenmatrices (JADE) can be employed.
Y=W·X (11)
The independent component matrix Y corresponds to an estimated value of the unknown component S. Accordingly, the expression (12) below can be obtained and the expression (13) below can be obtained by transforming the expression (12).
X=·Y (12)
Â=X·Y
+ (13)
In the expression, ̂A represents an estimated mixed matrix of A and Y+ indicates a pseudo inverse matrix of Y.
The estimated mixed matrix ̂A (this notation is made because of the restriction of characters used in the specification and actually means the character with a symbol on the left side of the expression (13), and the same applies to other characters) obtained by the expression (13) can be represented by the following expression.
In Step S110 of
After the process of Step S110 is finished, the CPU 10 performs a process of calculating the independent component matrix Y based on the separation matrix W and the spectrum data X for each sample obtained in the process 2 and stored in the hard disk drive 30 in advance (Step S120). The calculation process is a process of performing calculation according to the expression (11) above. In the processes of Steps S110 and 120, the CPU 10 functions as the independent component matrix calculation unit 432 of
Next, the CPU 10 performs a process of calculating the estimated mixed matrix ̂A based on the spectrum data X for each sample stored in the hard disk drive 30 in advance and the independent component matrix Y calculated in Step S120 (Step S130). The calculation process is a process of performing calculation according to the expression (13) above.
The estimated mixed matrix ̂A is obtained by the processes up to Step S130. That is, an element (estimated mixing coefficient) ̂aij of the estimated mixed matrix ̂A is obtained. That is, the estimated mixing coefficient ̂aij corresponds to the inner product value P calculated in (D) to (F) in
In Step S140, the CPU 10 acquires a correlation (degree of similarity) between target component contents C1, C2, . . . , Cn measured by the process 3 and components of respective columns included in the estimated mixed matrix ̂A (hereinafter, referred to as a mixing coefficient vector ̂α) calculated in Step S130. Specifically, a correlation between the target component content C (C1, C2, . . . , Cn) and the mixing coefficient vector ̂α1 in the first column (̂a11, ̂a21, . . . , ̂an1) is acquired, a correlation between the target component content C (C1, C2, . . . , Cn) and the mixing coefficient vector ̂α2 in the second column (̂a12, ̂a22, . . . , ̂an2) is acquired, and then a correlation with respect to the target component content C related to each column is sequentially acquired.
As the index indicating the degree of the correlation, a correlation coefficient R following the expression below can be used. The correlation coefficient R is referred to as Pearson's product-moment correlation coefficient.
−C and −̂αk each independently represent a target component amount and an average value of elements of a vector ̂αk.
As a result of Step S140 in
The selection in Step S150 is to select one column from among a plurality of columns when considered with the table TB of
The process 5 is a process of calculating a regression formula and performed using the computer 100 in the same manner as in the process 4. In the process 5, the computer 100 performs the process of calculating a regression formula of a calibration curve. Moreover, the process 5 may be performed after the data up to the process 4 is moved to another computer or a device.
C=u·P+v (16)
Here, C represents a target component content, P represents an inner product between measurement data and an independent component, and u and v represent an integer.
After the process in Step S210 is performed, the CPU 10 stores the constants u and v of the regression formula acquired in Step S210 and an independent component Yk corresponding to a target component order k (
The calibration method of a target component will be described below. The subject is set to be configured of the same components as those of a sample used when a calibration curve is generated. Specifically, the calibration method of a target component is performed using a computer. Moreover, the computer here may be the computer 100 used when a calibration curve is generated or another computer.
Xp={Xp1, Xp2, . . . , Xp3} (17)
In the process of Step S310, the CPU 10 functions as the subject observation data acquisition unit 510 of
After the process of Step S320 is performed, pre-processing is performed with respect to the observation data (absorbance spectrum Xp) of the subject obtained in Step S310 (Step S330). As the pre-processing, it is preferable to perform the same processing as that (that is, the normalization process performed by the first pre-processing unit 450 and the whitening process performed by the second pre-processing unit 460) used in the process 4 (more specifically, Step S110 of
Thereafter, the CPU 10 acquires the inner product value P between the independent component included in the calibration dataset DS2 and the pre-processed spectrum (pre-processed observation data) obtained in Step S330 (Step S340). The process of Step S340 corresponds to the process of (B) and (C) in
In the processes of Steps S330 and S340, the CPU 10 functions as the mixing coefficient calculation unit 530 of
Next, the CPU 10 acquires the content C of the target component by reading the constants u and v of the regression formula included in the calibration dataset DS2 from the hard disk drive 30 (the non-volatile storage device 550 in
Further, in the present embodiment, the content C acquired in Step S350 is set as the content of the target component of the subject, but, alternatively, the content C acquired in step S350 is corrected by a normalization coefficient used for normalization in Step S330 and then the corrected value may be used as the content to be acquired. Specifically, an absolute value (gram) of the content may be acquired by multiplying the content C by a standard deviation. According to the configuration, depending on the kind of target component, it is possible to make the content C have improved precision.
According to the above-described calibration method, the content of the target component can be acquired with high precision from one spectrum which is a measured value of the subject.
Hereinafter, various algorithms used in the first pre-processing unit 450, the second pre-processing unit 460, and the independent component analysis processing unit 470 shown in
As the first pre-processing performed by the first pre-processing unit 450, the noise addition process, standard normal variate transformation (SNV), and project on null space (PNS) can be used. The contents and the effects of the noise addition process have been described in
The SNV is given by the expression below.
Here, z represents processed data, x represents data to be processed (in the present example, the absorbance spectrum), xave represents an average value of data x to be processed, and σ represents a standard deviation of data x to be processed. As a result of the standard normal variate transformation, normalized data z in which the average value is 0 and the standard deviation is 1 is obtained.
When the PNS is performed, it is possible to decrease the base line fluctuation included in the data to be processed. In measurement of the data to be processed (the absorbance spectrum in the present embodiment), variation referred to as the base line fluctuation, for example, fluctuation in the average values of data for each piece of measurement data is generated between data due to various factors. Therefore, before the independent component analysis (ICA) is performed, it is preferable to remove the factors of fluctuation. The PNS can be used as pre-processing capable of decreasing the base line fluctuation in the data to be processed. Particularly, in regard to measurement data of an absorbance spectrum including an infrared region and a reflected light spectrum, since such base line fluctuation is large, application of the PNS is advantageous. Hereinafter, a principle in which the base line fluctuation included in data obtained by measurement (also simply referred to as “measurement data x”) is removed by the PNS will be described. Further, as a typical example, a case where the measurement data is an absorbance spectrum including an infrared region or is a reflected light spectrum will be described. In this case, in regard to other kinds of measurement data (for example, audio data or the like), the PNS can be applied in the same manner as described above.
In an ideal system, the measurement data x (data x to be processed) is represented by the expression below using m (m is an integer of 2 or greater) independent components si (i=1 to m) and respective mixing ratios ci.
Here, A represents a matrix (mixed matrix) formed by a mixing ratio ci.
In the independent component analysis (ICA), the process is performed on the assumption of this model. However, various fluctuation factors (the state of a sample or a change in the measurement environment) are present in actual measurement data. For this reason, a model inconsideration of the fluctuation factors, a model expressing the measurement data x using the expression below can be considered.
Here, the parameter b represents the amount of fluctuation in an amplitude direction of a spectrum, the parameter a represents the amount of constant base line fluctuation E (also referred to as “fluctuation in the average value”), the parameters b1, . . . , bg represent the amount of g (g is an integer of 1 or greater) fluctuations f1 (λ) to fg(λ) depending on the wavelength, and ε represents a fluctuation component other than those described above. Further, the constant base line fluctuation E is given by “E={1, 1, 1, . . . , 1}T (T on the right side represents transposition) and is a constant vector whose data length is equivalent to a data length N (the number of segments in the wavelength bandwidth) of the measurement data x. As a variable λ indicating a wavelength, N integers from 1 to N are used. That is, the variable λ corresponds to an ordinal number of the data length N (N is an integer of 2 or greater) of the measurement data x. At this time, fluctuations f1(λ), . . . , fg(λ) depending on the wavelength are given by “f1(λ)={f1(1), f1(2), . . . , f1(N)}T, . . . , fg(λ)={fg(1), fg(2), . . . , fg(N)}T.” Since these fluctuations are error factors in the ICA or calibration, it is desirable to remove the fluctuations in advance.
As the function f(λ), it is preferable to use one variable function in which the function f(λ) value monotonically increases according to an increase of λ in the range in which the value of λ is 1 to N. In the project on null space, the fluctuation included in the measurement data can be further reduced when a function other than an exponential function λα of λ in which an exponent α is an integer is used.
As a method of determining the functional types of a preferable function f(λ) and the number g thereof, experimental trial and errors may be employed or existing parameter estimation algorithms (for example, an expectation maximization method (EM) algorithm) may be used.
In the PNS, when a space formed of the above-described respective base line fluctuation component E and f1(λ) to fg(λ) is considered and the measurement data x is projected to a space (null space) without having these fluctuation components, data with reduced base line fluctuation component E and f1(λ) to fg(λ) can be obtained. As a specific operation, data z processed by the PNS is calculated by the expression below.
Here, P+ represents a pseudo inverse matrix of P. ki is obtained by projecting a constituent component si of the expression (20) to a null space without having fluctuation components. Further, ε* is obtained by projecting a fluctuation component ε of the expression (20) to a null space.
Moreover, when normalization (for example, the SNV) is performed after the process of the PNS, it is possible to remove the influence of a fluctuation amount b in the amplitude direction of the spectrum in the expression (20).
When the ICA is performed on data pre-processed by the PNS, the obtained independent component becomes an estimated value of a component ki of the expression (21) and becomes a value different from the true constituent component si. However, since the mixing ratio ci is not changed from the original value in the expression (20), the mixing ratio ci does not affect the calibration process ((A) to (D) in
Further, the details of the PNS are described in “Extracting Chemical Information from Spectral Data with Multiplicative Light Scattering Effects by Optical Path-length Estimation and Correction” written by Zeng-Ping Chen, Julian Morris, and Elaine Martin, 2006.
As the second pre-processing performed by the second pre-processing unit 460, the principal component analysis (PCA) and the factor analysis (FA) can be used.
In a general technique of the ICA, dimension compression of data to be processed or decorrelation is performed as the pre-processing. Since a transformation matrix to be acquired by the ICA is restricted to an orthogonal transformation matrix by the pre-processing, the calculation amount of the ICA can be reduced. Such pre-processing is referred to as “whitening,” and the PCA is used in many cases. The whitening using the PCA is described in Chapter 6 of “Independent Component Analysis” written by Aapo Hyvarinen, Juha Karhumen, and Erkki Oja, published by John Wiley & Sons, Inc., 2001 (“Independent Component Analysis,” in February 2005, published by publishing department of Tokyo Denkki University).
However, in the PCA, in a case where random noise is included in the data to be processed, an error may be generated in the measurement results due to the influence of the random noise. Here, in order to reduce the influence of the random noise, it is preferable to perform whitening using the factor analysis (FA) having robustness with respect to noise in place of the PCA. Hereinafter, the principle of whitening using the FA will be described.
As described above, generally in the ICA, a linear mixed model (the expression (19) above) which represents the data x to be processed as a linear sum of the constituent components si is assumed and the mixing ratio ci and the constituent components si are acquired. However, in actual data, random noise other than the constituent components si is added in many cases. Here, as a model in consideration with the random noise, a model expressing the measurement data x using the expression below can be considered.
x=A·s+ρ (22)
Here, ρ represents the random noise.
In addition, whitening in consideration with the noise mixed model is performed and then an estimate of the mixed matrix A and the independent component si can be obtained by performing the ICA.
In the FA of the present embodiment, it is assumed that the independent component si and the random noise ρ respectively follow normal distribution N(0, Im) and N(0, Σ). Moreover, as is generally known, the first parameter x1in normal distribution N (x1, x2) represents an expected value and the second parameter x2 therein represents a standard deviation. At this time, since the data x to be processed becomes a linear sum of variables following the normal distribution, the data x to be processed also follows the normal distribution. Here, when a covariance matrix of the data x to be processed is set as V[x], the normal distribution followed by the data x to be processed can be implemented as N(0, V[x]). At this time a likelihood function related to the covariance matrix V[x] of the data x to be processed can be calculated by the following procedures.
First, when it is assumed that the independent components si are orthogonal to each other, the covariance matrix V[x] of the data x to be processed is calculated by the expression below.
V[x]=E|xx
T
|=AA
T+Σ (23)
Here, Σ represents a covariance matrix of noise ρ.
In this manner, the covariance matrix V[x] can be represented by the mixed matrix A and the covariance matrix Σ of noise. At this time, a log-likelihood function L (A, Σ) is given by the expression below.
L(A, Σ)=−n/2{tr((AAT+Σ)−1C)+log(det(AA1+Σ))+m log 2π} (24)
Here, n represents the number of pieces of data x, m represents the number of independent components, the operator tr represents a trace of a matrix (sum of diagonal components), and the operator det represents a determinant. In addition, C represents a sample covariance matrix acquired by sample calculation from the data x and is calculated by the expression below.
The covariance matrix Σ of the mixed matrix A and the noise can be acquired by the maximum-likelihood method using the log-likelihood function L (A, Σ ) of the expression (24). A matrix which is hardly affected by the random noise ρ of the expression (22) is obtained as the mixed matrix A. This is the basic principle of the FA. Further, as the algorithm of the FA, there are various algorithms using an algorithm other than the maximum likelihood method. In the present embodiment, various kinds of FAs can be used.
The estimated value obtained by the FA is only a value of AAT and the influence of the random noise is reduced and the data can be decorrelated in a case where a mixed matrix A suitable for the value of AAT is determined, but a plurality of constituent components si cannot be uniquely determined because the degree of freedom of rotation remains. Meanwhile, the ICA is a process of reducing the degree of freedom of rotation of the plurality of constituent components si such that the plurality of constituent components si are orthogonal to each other. In the present embodiment, the values of the mixed matrix A acquired by the FA are used as a whitened matrix and arbitrary properties with respect to the remaining rotation are specified by the ICA. In this manner, after a robust whitening process is performed on the random noise, independent constituent components si which are orthogonal to each other can be determined by performing the ICA. In addition, as a result of such process, the calibration precision in regard to the constituent components si can be improved by reducing the influence of the random noise.
In the independent component analysis (ICA), as the index for separating independent components from each other, higher order statistics representing independence of separated data are used as the independence index. The kurtosis is a typical independence index. The ICA using the kurtosis as the independence index is described in Chapter 8 of “Independent Component Analysis” written by Aapo Hyvarinen, Juha Karhumen, and Erkki Oja, published by John Wiley & Sons, Inc., 2001 (“Independent Component Analysis,” in February 2005, published by publishing department of Tokyo Denkki University).
However, in a case where an outlier such as spike noise is present in the data to be processed, the statistics including the outlier are calculated as the independence index. For this reason, an error is generated between original statistics related to the data to be processed and calculated statistics and this causes degradation of separation precision. Here, it is preferable to use an independence index which is unlikely to be affected by the outlier in the data to be processed. It is possible to use the β divergence as the independence index having such characteristics. Hereinafter, the principle of β divergence as the independence index of the ICA will be described.
As described above, generally in the ICA, a linear mixed model (the expression (19) above) which represents the data x to be processed as a linear sum of the constituent components si is assumed and the mixing ratio ci and the constituent components si are acquired. An estimated value y of the constituent component s acquired by the ICA is represented by “y=W·y” using a separation matrix W. It is desired that the separation matrix W is an inverse matrix of the mixed matrix A.
Here, a log-likelihood function L (̂W) of an estimated value ̂W of the separation matrix W can be represented by the expression below.
Here, an element of an integration symbol Σ is a log likelihood in each data point x(t). The log-likelihood function L (̂W) can be used as an independence index in the ICA. The technique of β divergence is a method of transforming the log-likelihood function L (̂W) with respect to an outliner such as spike noise in data for the purpose of suppressing the influence of the outliner by acting an appropriate function on the log-likelihood function L (̂W).
In a case of using the β divergence as the independence index, first, the log-likelihood function L (̂W) is transformed by the expression below using a function Φβ selected in advance.
In addition, the function LΦ (̂W) is considered as a new likelihood function.
As the function Φβ for reducing the influence of an outliner such as spike noise, a function in which the value of the function Φβ is exponentially attenuated when the value of the log likelihood (value in parentheses of the function Φβ) becomes smaller is considered. As such a function Φβ, the following function can be used.
In this function, the function value with respect to each data point z (log likelihood in the expression (27) above) becomes smaller when the value of β becomes larger. The value of β can be empirically determined and can be set as, for example, approximately 0.1. Further, the function Φβ is not particularly limited to the function of the expression (28) and another function in which the function value with respect to each data point z becomes smaller when the value of β becomes larger can be used.
When such β divergence is used as the independence index, the influence of an outliner such as spike noise can be appropriately suppressed. When the likelihood function LΦ (̂W) as in the expression (27) is considered, a pseudo distance between probability distribution to be minimized in correspondence with the maximization of likelihood is the β divergence. When the ICA using such β divergence as the independence index is performed, the influence of an outliner such as spike noise is reduced and the calibration precision related to the constituent component si can be improved.
In addition, the ICA using the β divergence is described in “Robust Blind Source Separation by β-Divergence” written by Minami Mihoko and Shinto Eguchi, 2002.
The invention is not limited to examples or modification examples and can be implemented with various aspects within the range not departing from the scope of the invention, and the following modifications are possible.
In the embodiment, the element number m of the spectrum S of an unknown component is empirically and experimentally determined in advance, but the element number m of the spectrum S of an unknown component may be determined based on the information criterion known as the minimum description length (MDL) or the akaike information criteria (AIC). In a case of using the MDL or the like, the element number m of the spectrum S of an unknown component can be automatically determined by an operation from observation data of samples. In addition, for example, the MDL is described in “Independent component analysis for noisy data? MEG data analysis, 2000.”
In the embodiment, the subject as a target of the calibration process is configured of the sample components as those of a sample used when the calibration curve is generated, but unknown components other than the components which are the same as those of the sample used when the calibration curve is generated may be included in the subject. This is because the inner product between independent components is set as 0 and the inner product of independent components corresponding to the unknown component is considered as 0, the influence of the unknown component in a case of acquiring a mixing coefficient using the inner product can be ignored.
The computer used in the embodiment may be configured as a dedicated device. For example, the device shown in
In the embodiment, an input of a spectrum of the spectral reflectance related to a sample or a subject is performed by inputting the spectrum measured by a spectrometer, but the invention is not limited thereto. For example, a spectrum is estimated from a plurality of band images whose wavelength bandwidths are different from each other and the spectrum may be input. The band images can be obtained by imaging a sample or a subject using a multi-band camera including a filter capable of changing a transmission wavelength bandwidth.
Further, elements other than the elements described in the aspects from among the constituent elements in each of the examples and modification examples described above are additional elements and can be appropriately omitted.
Number | Date | Country | Kind |
---|---|---|---|
2014-207163 | Oct 2014 | JP | national |