WO2022270289A1 - 検体に含まれる成分の含有量比の推定方法、組成推定装置、及び、プログラム - Google Patents

検体に含まれる成分の含有量比の推定方法、組成推定装置、及び、プログラム Download PDF

Info

Publication number
WO2022270289A1
WO2022270289A1 PCT/JP2022/022813 JP2022022813W WO2022270289A1 WO 2022270289 A1 WO2022270289 A1 WO 2022270289A1 JP 2022022813 W JP2022022813 W JP 2022022813W WO 2022270289 A1 WO2022270289 A1 WO 2022270289A1
Authority
WO
WIPO (PCT)
Prior art keywords
matrix
sample
specimen
intensity distribution
components
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/JP2022/022813
Other languages
English (en)
French (fr)
Inventor
昌信 内藤
裕理 日比
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National Institute for Materials Science
Original Assignee
National Institute for Materials Science
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National Institute for Materials Science filed Critical National Institute for Materials Science
Priority to JP2023529793A priority Critical patent/JP7628731B2/ja
Priority to US18/569,560 priority patent/US20240297031A1/en
Priority to EP22828206.7A priority patent/EP4361624B1/en
Publication of WO2022270289A1 publication Critical patent/WO2022270289A1/ja
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • H—ELECTRICITY
    • H01—ELECTRIC ELEMENTS
    • H01J—ELECTRIC DISCHARGE TUBES OR DISCHARGE LAMPS
    • H01J49/00—Particle spectrometers or separator tubes
    • H01J49/0027—Methods for using particle spectrometers
    • H01J49/0036—Step by step routines describing the handling of the data generated during a measurement
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/20—Identification of molecular entities, parts thereof or of chemical compositions
    • H—ELECTRICITY
    • H01—ELECTRIC ELEMENTS
    • H01J—ELECTRIC DISCHARGE TUBES OR DISCHARGE LAMPS
    • H01J49/00—Particle spectrometers or separator tubes
    • H01J49/0009—Calibration of the apparatus

Definitions

  • the present invention relates to a method for estimating the content ratio of components contained in a sample, a composition estimation device, and a program.
  • a method for estimating the composition of a sample which is a mixed system of multiple components
  • a method of separating each component in the sample by some method and identifying and quantifying each fraction component is known. It is As one of the identification methods, a mass spectrum obtained by a mass spectrometer has been used.
  • Non-Patent Document 1 a method of reducing the dimensions of these multidimensional data by principal component analysis and searching for features has been used.
  • mass spectra in which the peak intensity depends on the ionization efficiency.
  • ionization efficiency which is an unknown variable for each compound.
  • isotope samples with known sample amounts and the same ionization efficiency as standard substances. there were.
  • Non-Patent Document 1 loses a lot of information of the original multidimensional data in the process of dimensionality reduction, and although it can perform qualitative analysis of components, it cannot perform quantitative analysis. I didn't.
  • the ratio of the magnitude (area) of the mass spectrum obtained is the ratio of the content of each component in the mixture. was not reflected, and this was one of the factors hindering quantitative analysis by mass spectroscopy.
  • ADI ambient desorption ionization
  • the present invention aims to provide a method for estimating the content ratio of components that can estimate the composition of the components of an unknown mixture even from a mass spectrum obtained under ambient conditions. Make it an issue.
  • Another object of the present invention is to provide a composition estimation device and a program.
  • [1] A method for estimating the content ratio of the component in an estimation target specimen containing at least one component selected from K types of components, where K is an integer of 1 or more, wherein the K Preparing K or more training specimens having different compositions that contain at least one component selected from the following types of components, and a background specimen that does not contain the above components, the estimation target specimen, and the learning While heating each sample contained in the sample set including the sample for use and the background sample, the gas components generated by thermal desorption and/or thermal decomposition are sequentially ionized, and the mass spectrum is continuously observed.
  • Extracting the noise component contained in the intensity distribution matrix by, correcting the intensity distribution matrix so as to reduce the influence of the noise component, obtaining a corrected intensity distribution matrix, and converting the corrected intensity distribution matrix to the specimen Divide into small matrices corresponding to each, express each of the specimens in the vector space using the small matrices as feature vectors, and set a K-1-dimensional simplex containing all of the feature vectors, Determining K endmembers in a K-1-dimensional simplex, calculating Euclidean distances between the K endmembers and the feature vector of the estimated target specimen, respectively, and calculating the Euclidean distance ratio, and estimating the content ratio of the component in the estimation target specimen, and when the K is 3 or more, the learning for A method wherein at least one of the specimen feature vectors is located or wherein the training specimen comprises at least one of the endmembers.
  • the determination of the endmembers performs a non-negative matrix factorization of the corrected intensity distribution matrix to obtain a matrix representing the weight fractions of the K types of the components in the specimen; A method according to [6], by a second NMF process that decomposes into products with a matrix representing the fragment presence intensities of .
  • a method according to [6] by a second NMF process that decomposes into products with a matrix representing the fragment presence intensities of .
  • the method of [6] or [7] wherein the training sample does not contain the endmember.
  • obtaining the corrected intensity distribution matrix further includes performing intensity correction on the intensity distribution matrix.
  • the sample set further includes a calibration sample, the calibration sample includes all of the K types of components, and has a known composition, and the intensity correction is performed to convert the intensity distribution matrix into The method of [9], comprising normalizing.
  • the intensity correction includes allocating at least a portion of the intensity distribution matrix by the product of the mass of the corresponding specimen in the specimen set and an environmental variable; The method according to [9], wherein the environmental variable is a variable representing the effect of the component on ionization efficiency during the observation.
  • the environmental variable is a compound with a molecular weight of 50 to 1500 contained in the atmosphere at the time of the observation, or an organic low molecular weight compound with a molecular weight of 50 to 500 contained in each of the specimen sets in a predetermined amount.
  • a composition estimation device for estimating the content ratio of the component in an estimation target specimen containing at least one component selected from K types of components, where K is an integer of 1 or more, Included in a sample set containing K or more training samples having different compositions, background samples not containing the above components, and the estimation target sample, which contain at least one component selected from the above K types of components
  • a mass spectrometer that successively ionizes gas components generated by thermal desorption and/or thermal decomposition while heating each specimen, and continuously observes the mass spectrum, and information for processing the observed mass spectrum.
  • the information processing device stores the mass spectrum obtained for each heating temperature in each row, obtains the two-dimensional mass spectrum for each specimen, and converts at least two of the two-dimensional mass spectra
  • a data matrix creation unit that collectively converts the above into a data matrix
  • NMF processing that performs non-negative matrix factorization of the data matrix and decomposes into the product of a normalized basis spectrum matrix and its intensity distribution matrix
  • NMF A processing unit extracts a noise component included in the intensity distribution matrix by canonical correlation analysis between the base spectral matrix and the data matrix, corrects the intensity distribution matrix so as to reduce the influence of the noise component, and corrects the a correction processing unit that creates an intensity distribution matrix; divides the corrected intensity distribution matrix into submatrices corresponding to each of the specimen sets; a vector processing unit, an end member determination unit for setting a K-1 dimensional simplex containing all of the feature vectors, and determining K end members in the K-1 dimensional simplex; and the feature vector of the estimation target specimen,
  • the sample for use comprises at least one of said endmembers.
  • the end member determination unit determines the end members based on discriminant labels attached to the learning samples.
  • the end member determination unit determines the end members by an algorithm that determines vertices so that the K-1 dimensional simplex has the minimum volume.
  • the end member determination unit performs non-negative matrix factorization on the corrected intensity distribution matrix, and obtains a matrix representing the weight fraction of the K types of the components in the specimen, and the K types of individual components.
  • the specimen set further includes a calibration specimen, the calibration specimen contains all of the K types of components, and has a known composition, and the intensity correction is performed to convert the intensity distribution matrix to The composition estimation device according to [17], including normalizing.
  • the intensity correction includes assigning at least a portion of the intensity distribution matrix by the product of the mass of the corresponding analyte in the analyte set and an environmental variable;
  • a mass spectrometer that sequentially ionizes gas components generated by thermal desorption and/or thermal decomposition while heating each sample contained in a sample set and continuously observes the mass spectrum, and the observed mass spectrum.
  • a data matrix creation function that converts at least two or more of the above into a data matrix collectively, and NMF processing that decomposes the data matrix into a product of a normalized basis spectrum matrix and its intensity distribution matrix by non-negative matrix factorization extracting the noise component contained in the intensity distribution matrix by the NMF processing function and the canonical correlation analysis between the base spectrum matrix and the data matrix, and correcting the intensity distribution matrix so as to reduce the influence of the noise component.
  • a correction processing function that creates a post-correction intensity distribution matrix by dividing the post-correction intensity distribution matrix into submatrices corresponding to each of the specimens, and using the submatrices as feature vectors, the specimens in the vector space and an end member determination function for setting a K-1 dimensional simplex containing all of the feature vectors and determining K endmembers in the K-1 dimensional simplex, Euclidean distances between the K end members and the feature vector of the estimation target specimen are calculated, and the content ratio of the component in the estimation target specimen is estimated from the ratio of the Euclidean distances.
  • the K is 3 or more, at least one of the feature vectors of the learning specimen is located in each of the regions outside the hypersphere inscribed in the K-1-dimensional simplex, or , the program, wherein the training sample comprises at least one of the endmembers.
  • the method of the present invention heats each sample contained in a sample set, sequentially ionizes gas components generated by thermal desorption and / or thermal decomposition, and obtains a mass spectrum. Observing continuously, obtaining a two-dimensional mass spectrum for each specimen in which the obtained mass spectrum for each heating temperature is stored in each row, and normalizing the data matrix summarizing these by non-negative matrix factorization (NMF)
  • NMF non-negative matrix factorization
  • the noise component contained in the intensity distribution matrix is extracted and corrected by canonical correlation analysis, and the corrected intensity distribution matrix containing the feature vector of each specimen is obtained.
  • One of the characteristic points is the subsequent quantitative calculation.
  • the present method determines the mixture ratio (composition ) can be estimated.
  • the composition of the sample can be estimated more accurately and easily.
  • the endmembers can be determined from the label information, so the composition of the sample can be estimated more easily.
  • the endmembers can be determined by vertex component analysis of the feature vector of the training specimen, and the composition of the specimen can be estimated. can.
  • the composition of the sample can be determined by determining the end member by using an algorithm that determines the vertex so that the K-1 dimensional simplex has the minimum volume. can be estimated. It is excellent in that it does not require an end member as a training sample.
  • the post-correction intensity distribution matrix is the post-intensity correction of the intensity distribution matrix, the accuracy of the estimation is further enhanced.
  • the composition of the sample can be estimated with higher accuracy using a simpler procedure.
  • the composition of the specimen can be estimated with high accuracy.
  • a compound with a molecular weight of 50 to 1500 contained in a predetermined amount in the atmosphere at the time of observation, or an organic low molecule with a molecular weight of 50 to 500 contained in a predetermined amount in each sample of the sample set If it is the total value of the peaks of the mass spectrum of the compound, even if the sample set does not include the calibration sample, a simpler operation will yield a more accurate composition estimation result.
  • FIG. 4 is a flow chart of a composition estimation method according to an embodiment of the present invention
  • FIG. 10 is an image diagram showing a K ⁇ 1 dimensional simplex (two-dimensional simplex, triangle) when K is 3
  • 1 is a hardware configuration diagram of a composition estimation device according to an embodiment of the present invention
  • FIG. 1 is a functional block diagram of a composition estimating device according to an embodiment of the present invention
  • FIG. 10 is a diagram showing the original spectrum (A) and the reconstructed spectrum (B) after processing for "MSE_calib"; Part C ⁇ (MSE_Calib) of the C matrix after NMF decomposition and S T after removing background and impurity spectra by CCA-Filter.
  • FIG. 10 is a diagram showing the original spectrum (upper, raw of MSE_calib) for "MSE_calib" and the reconstructed spectrum after performing the above processing with the CCA-filter removed. It is a composition estimation result when correction by an environmental variable is not performed.
  • FIG. 10 is an image diagram showing a K ⁇ 1 dimensional simplex (two-dimensional simplex, triangle) when K is 3; It is an explanatory diagram of the mechanism of the end member determination method in the second embodiment of the present invention.
  • FIG. 11 is a diagram showing estimation results of Example 2;
  • FIG. 11 is a diagram showing estimation results of Example 3;
  • component means a substance that constitutes each sample, and estimation target samples, learning samples, and calibration samples mean mixtures made up of combinations of components.
  • Presume target sample is a sample whose composition should be estimated, and both the number of components contained and the content of each component may be unknown.
  • the "learning sample” means a sample necessary for determining K endmembers of the K-1 dimension simple substance described later, and these samples include at least one component selected from K components. contained, each with a different composition.
  • the number of components contained is not particularly limited, and may be from 1 to K types.
  • the learning sample may contain the same composition as the estimation target sample. That is, the estimation target sample and the learning sample may be the same, but the learning samples are different.
  • End member means a vector corresponding to the vertex of a K-1 dimensional simplex, and this end member is a sample (learning sample or virtual sample) consisting of only one component of K types. Corresponds to the feature vector.
  • a “background sample” means a sample that is different from the estimation target sample and the learning sample and that does not contain any of the above components.
  • a “calibration sample” means a sample that is different from the estimation target sample, the learning sample, and the background sample, contains all K types of components, and has a clear composition.
  • FIG. 1 is a flow chart of a composition estimation method (hereinafter also referred to as “this method”) according to an embodiment of the present invention.
  • an estimation target specimen, K or more learning specimens with mutually different compositions, and background specimens are prepared (step S101).
  • the sample set may further include a calibration sample.
  • the estimated target sample contains at least one of the K types of components. If the analyte of interest contains one type of component, the purity of that component can be estimated by this method. In this case, the presumed target specimen includes the above-mentioned "ingredient" and other impurities.
  • the putative target analyte is a mixture of components.
  • the upper limit is not particularly limited.
  • K is preferably 1000 or less.
  • the number of learning samples is related to the number of components K, and must be K or more. If it is K or more, the number of learning samples is not particularly limited, but if it is more, the accuracy of determining the vertices of the K-1-dimensional simple substance described later tends to be higher, and as a result, more accurate composition estimation. Easy to get results.
  • the ratio of the number of learning samples to the number of components K is 1 or more, preferably 2 or more, more preferably 3 or more, further preferably 5 or more, and 10 or more.
  • the upper limit is not particularly limited, but generally 1000 or less is preferable.
  • the learning samples differ in composition from each other.
  • the composition is determined by the type and number of components contained and the content of each component.
  • the same composition means that the type, number, and content of the contained components are the same, and if any one of them is different, it will be in a state of "different compositions.”
  • the learning samples have different compositions
  • other information that is, the types of components contained, the number of components, and the contents of each component, is used for estimating the composition of the estimation target sample.
  • the composition of the training sample may be unknown.
  • the difference between the estimation target sample and the learning sample is that the learning sample must have different compositions.
  • the learning samples may include those having the same composition as the estimation target sample.
  • the training sample may contain end members. That is, the learning sample may contain one type of component selected from K types of components (simple, pure). When the training sample contains end members, it is possible to more easily and accurately estimate the composition.
  • the end-member learning sample may include a discriminant label. It is preferable that the training sample contains a discriminant label, because it makes it easier to determine the endmembers.
  • each of the estimation target specimens is heated to sequentially ionize gases generated by thermal desorption and/or thermal decomposition to acquire a mass spectrum (step S102).
  • the method of obtaining a mass spectrum is not particularly limited, a method of performing mass spectrometry on a sample under ambient conditions without pretreatment is preferred.
  • a mass spectrometer called DART-MS, which combines an ion source called DART (Direct Analysis in Real Time) ion source and a mass spectrometer, is known.
  • the mass spectrometer is not particularly limited, and is preferably capable of accurate mass spectrometry, and may be either a quadrupole type, a time-of-flight (TOF) type, or the like.
  • each sample contained in the sample set is sequentially heated at a heating rate of 50 ° C./min, and heated to 50 to 550
  • the mass spectrum obtained for each heating temperature is stored in each row to obtain a two-dimensional mass spectrum for each specimen, and at least two or more of these two-dimensional mass spectra are collectively converted into a data matrix (step S103 ).
  • step S102 mass spectra are continuously acquired at predetermined temperature rise intervals. These mass spectra may be used to create a data matrix as they are, or may be averaged for each predetermined heating temperature range. The amount of data can be compressed by averaging the mass spectra for each predetermined heating temperature range into one line. Such a temperature rise range is, for example, about 10 to 30.degree.
  • the peak intensity may be normalized.
  • a method of normalization for example, a method of normalizing so that the sum of squares of peak intensities becomes 1 can be used.
  • the number of two-dimensional mass spectra used to create the data matrix X is not particularly limited as long as it is two or more, it is preferable to use the two-dimensional mass spectra of all samples included in the sample set.
  • part or all of the two-dimensional mass spectrum obtained by the measurement two or more times may be used to create the data matrix X. .
  • NMF non-negative matrix factorization
  • the spectral data matrix is generally
  • the data matrix X is an N ⁇ M data matrix formed by stacking N spectra with M channels
  • S is a base spectrum matrix
  • C is its intensity distribution matrix. If the number of all elements included in the data set is K1, then C is an N ⁇ K1 matrix and S is an M ⁇ K1 matrix. Since the actual measurement contains noise,
  • the NMF processing in this step is performed by alternately updating the column vectors of C and S according to Hierarchical alternating least square (HALS), which is easy to incorporate the constraints described below. is preferably solved as a convex optimization problem for each column of .
  • HALS Hierarchical alternating least square
  • a mass spectrometer generally has a resolution of about 0.1 to 0.001 m/z. Therefore, situations in which the m/z of different gas components match to the nearest decimal point and occupy the same channel are rare. Therefore, it is reasonable to assume some degree of orthogonality between the column vectors of the basis spectral matrix S.
  • Threshold is a hyperparameter given between 0.9 and 0.99. It is preferable to merge similar S in a similar manner, and also to superimpose the corresponding column vectors of S according to the weight of C, even when multiple column vectors of C are parallel. i.e.
  • the element ⁇ k of is a parameter following an inverse gamma distribution representing the importance of the component k in the data. Since we assumed a Gaussian distribution for the noise
  • a is a hyperparameter such that a-1 is a very small positive number
  • b is the average value ⁇ X of the data matrix.
  • ⁇ k is the Lagrangian undetermined multiplier
  • the soft orthogonality condition can freely adjust the degree of orthogonality between complete orthogonality and no orthogonality constraint by multiplying ⁇ k by a hyperparameter w of [0, 1].
  • the final update formula is
  • the base spectrum matrix and the data matrix are subjected to canonical correlation analysis to obtain a corrected intensity distribution matrix (step S105).
  • the noise component contained in the intensity distribution matrix is extracted by canonical correlation analysis between the base spectrum matrix and the data matrix, the intensity distribution matrix is corrected to reduce the influence of the noise component, and the corrected intensity distribution matrix is get
  • NMF is a low-rank approximation of the data matrix, even if the component k does not actually exist in the i-th spectrum, if the approximation is better in the sense of least squares, then C ik >0 and In many cases such C ik is very small and often does not matter in NMF analysis.
  • C jk >0 which is present in trace amounts in the j-th spectrum
  • C ik >0 which is an NMF artifact
  • substitute 0 for C ik which is a ghost peak
  • Canonical Correlation Analysis is used to solve the above.
  • This step is one of the features of this method, and the present inventors named this step Canonical Correlation Analysis (CCA) Filter.
  • CCA Canonical Correlation Analysis
  • the CCA filter scans sample-wise whether each component of the base spectrum S T output from NMF was actually included in the original data, and if no peak pattern similar to the original data is found, By erasing from the spectrum of that sample. If the spectrum generated from the L-th sample is stored in rows 1 to 1' of the data matrix X, the submatrix extracting that portion is written as X ⁇ (L) . , l'-2, l'-1 ⁇ .
  • CCA is conventionally known as a technique for examining the correlation between data of two systems, and is formulated as follows.
  • a matrix of N p-dimensional vectors is
  • a matrix consisting of N q-dimensional vector data is
  • the CCA problem setting is to project Y and Z onto a certain axis and find the solution (a T ,b T ) T that maximizes the correlation coefficient ⁇ .
  • the noise that can be removed by the CCA filter includes, for example, (1) NMF artifacts (so-called false peaks), (2) peaks derived from trace chemical substances that slip into the surrounding environment of the device (background), (3) However, there are peaks (contamination) derived from substances that may distort the nature of the analysis if included in the analysis.
  • Contaminants include, but are not limited to, solvents used in sample preparation and loading.
  • This step may further include a step of intensity correction of the intensity distribution matrix.
  • DART-MS primarily ionizes water and ammonia in the air with helium ions, which secondarily ionizes the sample. Therefore, the peak intensity may change depending on atmospheric conditions such as humidity. Also, the sample mass to be loaded is usually different for each sample, which also affects the peak intensity. Also, the ionization efficiency of each component contained in the sample is generally different, which also affects the peak intensity. By intensity correction of the intensity distribution matrix, the above effects can be reduced and more accurate quantification results can be obtained.
  • the intensity correction can be performed by normalizing the intensity distribution matrix.
  • Normalizing the intensity distribution matrix means normalizing the C matrix for each specimen, and more specifically means an operation of adding all elements to 1 for each specimen. If we divide the C matrix sample-wise and denote the C matrix derived from sample L as C ⁇ (L) ,
  • the ionization efficiency differs for each compound, and the thermal decomposition efficiency also differs, and the amount of compound and the peak intensity are not proportional.
  • the peak intensity ratio including the ionization efficiency of component k by NMF is already reflected in S :kT , and the intensity ratio between spectra is reflected in C :k . Thanks to this hierarchical data structure formed by NMF, more accurate quantitative compositional analysis is possible without obtaining unknown ionization efficiencies, as long as peak intensity ratios at the endmembers are taken into consideration.
  • the specimen set does not include calibration specimens
  • more accurate quantitative analysis is possible by assigning at least part (preferably all) of the intensity distribution matrix by the product of the mass of the corresponding specimen and the environmental variable.
  • the environmental variable is a variable representing the influence of each component on the ionization efficiency during observation.
  • a component with a high ionization efficiency has a high peak intensity
  • a component with a low ionization efficiency has a low peak intensity.
  • the peak intensity generally increases as the mass increases, and the peak intensity decreases as the mass decreases.
  • the C matrix is corrected as follows so that the intensity information between samples held by the intensity distribution matrix can be used without normalizing the intensity distribution matrix.
  • m L is the mass of sample L
  • I L is the peak intensity of (poly)dimethylsiloxane generated from the oil pump attached to the DART-MS, observed in the entire temperature range at 610 m/z.
  • IL reflects the standby state at the time of measurement, and there is a difference of about three times in peak intensity between rainy days and sunny days. corrected by this
  • I L is 610 m/z derived from (poly)dimethylsiloxane, but is not limited to the above, and is derived from a compound with a molecular weight of 50 to 1500 contained in a predetermined amount in the atmosphere at the time of observation. and/or derived from an organic low-molecular-weight compound having a molecular weight of 50 to 500 contained in a predetermined amount in each of all samples.
  • the corrected intensity distribution matrix is divided into submatrices corresponding to the respective specimens included in the specimen set, and each specimen is expressed in the vector space as a feature vector (step S106 ). Furthermore, a K-1 dimensional simplex containing all of these feature vectors (each corresponding to each specimen) is set, and K endmembers in the K-1 dimensional simplex are determined (step S107).
  • the training specimen contains an end member and the training specimen has a discriminative label (a label for distinguishing that the training specimen is an "end member"), based on that label End members are determined.
  • the training sample contains endmembers and the training sample is not labeled discriminatively, i.e. the training sample contains endmembers, but which training sample is the endmember If is unknown, the endmembers are determined by Vertex Component Analysis of the feature vector.
  • the endmembers are determined by an algorithm that determines the vertices so that the K-1 dimensional simplex has the minimum volume. It should be noted that the fact that the K ⁇ 1 dimensional simplex has the minimum volume means the minimum area in the case of a two-dimensional simplex (triangle).
  • At least one of the feature vectors of the training specimen is located in each of the regions outside the hypersphere inscribed in the K-1-dimensional simplex. (when K is 3 or more), the endmembers can be determined with high accuracy.
  • FIG. 2 is an image diagram showing a K-1 dimensional simplex (2-dimensional simplex, triangle) when K is 3.
  • the K-1 dimensional simplex 10 is a triangle whose vertices are the end members 13, 14 and 15 determined by the training specimens 16, 17 and 18.
  • the algorithm for determining the vertices so that the K-1 dimensional simplex has the minimum volume is not particularly limited, and a known algorithm can be applied.
  • Examples of such algorithms include the Minimum Volume Constraint-NMF method, the Minimum Distance Constraint-NMF method, the Simplex identification via split augmented Lagrangian method, and the Simplex volume minimization method.
  • the Euclidean distance between the K endmembers and the feature vector of the sample to be estimated is calculated, and the content ratio of each component in the sample is estimated (step S108).
  • the feature vector of the other training sample is K ⁇ It may be located in a region inside a hypersphere inscribed in a one-dimensional simplex.
  • FIG. 13 is an image diagram showing a K-1 dimensional simplex (two-dimensional simplex, triangle) when K is 3, as in FIG.
  • a K-1 dimensional simplex 60 is a triangle whose vertices are end members 64 and 65 determined by a learning sample 61 as an end member, another learning sample 62 and a learning sample 63 .
  • the difference from FIG. 2 is that the learning specimens 62 and 63 are located in the region inside the hypersphere inscribed in the K-1 dimensional simplex.
  • the other learning samples may be located on a hypersphere inscribed in the K-1 dimensional simple substance, or may be located in the outer region. may be located. That is, in this case, the positions of other learning samples may be arbitrary.
  • the other training sample contains 20% by mass or more of an endmember component (another endmember component) different from the endmember learning sample. It is preferably contained, and more preferably 40% by mass or more.
  • step S107 in the already described first embodiment if the learning sample does not contain an end member, that is, there is no learning sample located at the K-1 dimensional simplex vertex Contains how the end members (vertices) of a case are determined.
  • One of the features of the method according to this embodiment is that the determination of the endmembers is performed by non-negative matrix factorization of the matrix representing the fragment presence intensity for each sample, representing the weight fraction of the K types of components in the sample.
  • the second NMF process decomposes into the product of the matrix and the matrix representing the fragment presence intensity for each single unit of the K kinds of components.
  • FIG. 14 is a schematic diagram for explaining the mechanism.
  • Each sample is a mixture of K types of components (hereinafter, described as a mixture of "polymers"), but what can actually be observed is the spectrum of M types of fragment gases generated by thermal decomposition of each polymer. is a suitable linear summation of S.
  • the feature points according to the present embodiment have the coefficient matrix as given that, , and if it were possible to measure K polymer units, it would be If it is observed that and the idea that composition C can be obtained by performing a second NMF treatment.
  • step S104 (hereinafter also referred to as “first NMF processing”) and the NMF processing performed in determining the end members in step S107 (hereinafter “second NMF processing”) will be described. (also referred to as “processing”) will be described in detail. Note that, in principle, steps that are not described below are the same as the steps in the first embodiment.
  • N represents a non-negative or real N ⁇ M matrix.
  • a queue against represent the n-th row vector, m-th column vector, and (n, m)-element, respectively.
  • X T is the transposed matrix of X, represents the Frobenius norm.
  • Tr(X) represents a trace and diag(X) represents a diagonal matrix consisting of diagonal elements.
  • 1 N and 11 N represent N-dimensional vectors or (N,N) matrices with all elements being 1's.
  • I N represents an N-dimensional identity matrix.
  • N the number of samples
  • T the number of temperature band divisions
  • D the number of channels
  • M the number of base spectra
  • the first NMF was developed to be well adapted to MS (Mass Spectrometry) data interpretation by adding the following three points to the ARD-SO-NMF proposed by Shiga et al.
  • the overall probabilistic model is can be written as where p(S) is and p( ⁇
  • a,b) is an inverse gamma distribution parameterized by (a,b), i.e. and a is a hyperparameter that adjusts sparsity, a 1+10 ⁇ 16 , b is empirically estimated from the expected value E(A im ) of A im , be associated with This is further expressed as E(X id ) and and approximating E(X id ) by the average ⁇ X of X, From b is is determined. So the negative log-likelihood function is can be written as Since this is a downward convex function with respect to ⁇ , i.e. , the update formula for ⁇ is obtained. Up to this point, the derivation of ARD-SO-NMF reported by Shiga et al. was completely followed.
  • HALS hierarchical alternating least square
  • the objective function to be minimized is and the gradient with respect to a m and s m is and if these are set to zero, the update formula is becomes. Since s m is normalized for each update, we omitted the constant factor.
  • the non-negative constraint was satisfied by projecting into the non-negative quadrant at each update. Specifically, for example: can be projected by here means a vector with absolute values for each element.
  • the first step consists of the M-spectrum Let the obtained S be a spectral set similar to S m: and dissimilar spectral sets is to divide into This split is is fulfilled.
  • each spectral pair is as well as can be written as The similarity is evaluated by the correlation coefficient ⁇ , here, CCA's problem setting is can be written as the solution is given as a solution to the generalized eigenvalue problem, here is the eigenvector are arranged as column vectors in descending order of their eigenvalues. each eigenvalue is the correlation coefficient of y and z produced by linear combination with the corresponding (u * , v * ) as coefficients.
  • the background component was removed from the system by deleting the corresponding column vector of A and row vector of S.
  • the M'-component derived from the background has been removed (Hereinafter, MM' is replaced with M for simplicity, ) corrects the assigned intensity for the sample amount and the internal standard peak.
  • the submatrix part A (n) of A for sample n was converted into a one-dimensional vector may be used as the input for the second NMF as the feature vector of sample n, but since the temperature distribution of M-fragments is unnecessary information in composition analysis, all temperature bands are summed for each sample to obtain M - represent the fragment abundance (FA) may be used as an input. here for simplicity We will do a second NMF on .
  • NLS Non-negative least square
  • the volume term of the simplex spanned by the row vector of Let the non-orthogonality terms between the row vectors of Assuming ⁇ >0 and ⁇ [0,1] as weights, by is sought.
  • test data is projected onto the hyperplane of S and B.
  • S-hyperplane is obtained. i.e. Taking the sum over the temperature range, , by projection onto the B-hyperplane and normalization, is obtained. i.e.
  • the dataset and hyperparameters are as follows.
  • FIG. 3 is a hardware configuration diagram of the composition estimation device according to the embodiment of the present invention.
  • the composition estimation device 30 includes a mass spectrometer 31 and an information processing device 33 , and the mass spectrometer 31 includes a decompression pump 32 .
  • the information processing device 33 has a processor 34, a storage device 35, a display device 36, and an input device 37, and the mass spectrometer 31 and the information processing device 33 are configured to be able to exchange data with each other. .
  • the processor 34 is, for example, a microprocessor, processor core, multiprocessor, ASIC (application-specific integrated circuit), FPGA (field programmable gate array), and GPGPU (general-purpose computing on graphics processing units).
  • ASIC application-specific integrated circuit
  • FPGA field programmable gate array
  • GPGPU general-purpose computing on graphics processing units
  • the storage device 35 has a function of temporarily and/or non-temporarily storing various programs and data, and provides a working area for the processor 34 .
  • the storage device 35 is, for example, ROM (Read Only Memory), RAM (Random Access Memory), HDD (Hard Disk Drive), flash memory, SSD (Solid State Drive), and the like.
  • the input device 37 can receive input of various types of information and input of instructions to the composition estimation device 30 .
  • the input device 37 may be a keyboard, mouse, scanner, touch panel, or the like.
  • the display device 36 can display the status of the composition estimation device 30, the progress of the analysis, the composition estimation results, and the like.
  • the display device 36 may be a liquid crystal display, an organic EL (Electro Luminescence) display, or the like.
  • the display device 36 may be configured integrally with the input device 37 .
  • the display device 36 is a touch panel display, and the form which provides GUI(Graphical User Interface) is mentioned.
  • the information processing device 33 which includes a processor 34, a storage device 35, an input device 37, and a display device 36 that can communicate data with each other via a data bus, is typically a computer.
  • the mass spectrometer 31 is equipped with a vacuum pump 32 and is typically a mass spectrometer equipped with a DART ion source, a specimen heating device, and a time-of-flight mass spectrometer. Both the ion source and the mass spectrometer included in the mass spectrometer 31 are non-limiting examples, and the configuration of the mass spectrometer included in the present composition estimation apparatus is not limited to the above.
  • FIG. 4 is a functional block diagram of the composition estimation device 30 according to the embodiment of the present invention.
  • the composition estimating apparatus 30 has a mass spectrometer 41 that performs mass spectrometry on all specimens, and an information processor 42 that processes the mass spectrum obtained by the mass spectrometer 41 .
  • the mass spectrometer 41 is configured including the mass spectrometer 31, and is a function realized by the mass spectrometer 31 being controlled by the control unit 43, which will be described later.
  • the mass spectrometer 41 performs thermal desorption and/or thermal decomposition while heating the loaded specimen set (the specimen, the learning specimen, and the background specimen, symbol a in FIG. 4) to determine the gas components is ionized, mass spectroscopy is performed sequentially, and a mass spectrum is output (symbol c in FIG. 4).
  • an environmental factor caused by the environment in which the mass spectrometer 31 is installed (symbol b in FIG. 4, for example, (poly)dimethylsiloxane generated from the decompression pump 32) is also loaded into the mass spectrometer 31 at the same time. Affects the spectrum.
  • the control unit 43 is configured including the processor 34 .
  • the control unit 43 implements each function of the composition estimation device 30 by controlling each unit.
  • the storage unit 44 includes a storage device 35 .
  • the storage unit 44 stores in advance a program for controlling each unit, and provides a temporary or non-temporary storage area for data input/output by each unit.
  • the storage unit 44 is controlled by the control unit 43 and can read and write data.
  • the input unit 45 includes an input device 37. Also, the output unit 46 is configured including the display device 36 . By controlling these, the control unit 43 receives input from the user of the composition estimation device 30, stores it in the storage unit 44, transmits information stored in the storage unit 44 to each unit, It can be displayed to the user of the composition estimation device 30 via the output unit 46 .
  • the data matrix creation unit 47 is a function realized by executing the program stored in the storage unit 44 by the control unit 43 .
  • the data matrix creation unit 47 creates a data matrix X from the two-dimensional mass spectrum in which the mass spectrum (symbol c in FIG. 4) acquired for each heating temperature by the mass analysis unit 41 is stored in each row, and performs NMF processing to be described later. It is passed to the unit 48 (symbol d in FIG. 4).
  • the details of the data matrix X have already been explained.
  • the NMF processing unit 48 is a function realized by executing the program stored in the storage unit 44 by the control unit 43 .
  • the result is passed to the correction processing unit 49 (symbol e in FIG. 4).
  • NMF non-negative matrix factorization
  • the correction processing unit 49 is a function realized by executing the program stored in the storage unit 44 by the control unit 43 .
  • the correction processing unit 49 has a function of correcting the intensity distribution matrix by canonical correlation analysis between the base spectrum matrix and the data matrix created by the NMF processing unit 48 to create a post-correction intensity distribution matrix.
  • the created intensity distribution matrix after correction is passed to the vector processing unit 50 (symbol f in FIG. 4).
  • the correction processing unit 49 may further perform the already-described intensity correction in addition to the above-described correction by the canonical correlation analysis.
  • a calibration sample composition calibr sample composition
  • it is passed to the correction processing unit 49 via the input unit 45 (symbol k in FIG. 4).
  • CCA filter canonical correlation analysis
  • the vector processing unit 50 is a function realized by executing a program stored in the storage unit 44 by the control unit 43 .
  • the vector processing unit 50 divides the corrected intensity distribution matrix created by the correction processing unit 49 into sub-matrices corresponding to each of all the specimens, and expresses all the specimens in the vector space using these as feature vectors. This result is passed to the end member determination unit 51 (symbol g in FIG. 4).
  • the end member determination unit 51 is a function realized by executing a program stored in the storage unit 44 by the control unit 43 .
  • the end member determining unit 51 sets a K-1 dimensional simplex containing all the feature vectors generated by the vector processing unit 50, and determines K end members in the K-1 dimensional simplex. This result is passed to the content ratio calculator 52 (symbol h in FIG. 4).
  • the learning sample contains an end member and has a discriminant label
  • information on the discriminant label is passed from the input unit 45 to the end member determination unit 51 (symbol j in FIG. 4).
  • the discrimination label is not essential for end member determination (processing by the end member determination unit 51).
  • the content ratio calculation unit 52 is a function realized by executing the program stored in the storage unit 44 by the control unit 43 .
  • the content ratio calculation unit 52 calculates Euclidean distances between the K endmembers determined by the endmember determination unit 51 and the feature vectors of the estimation target specimen created by the vector processing unit 50, respectively. From the ratio of the Euclidean distance between the endmember and the feature vector of the specimen, the content ratio of K components in the specimen to be estimated is estimated.
  • the content ratio estimation result (indicated by symbol i in FIG. 4) is output from the output unit 46 . At this time, the result of estimating the composition of the learning sample may also be output.
  • composition estimation apparatus even solid samples that are difficult to dissolve in solvents and difficult to perform conventional composition analysis can be subjected to analysis without special pretreatment, and accurate quantitative analysis results can be obtained. be done.
  • This composition estimating apparatus can significantly reduce the time required to obtain results, particularly when applied to quality control of a sample that is a mixture of multiple components, investigation of causes of defects, and the like.
  • composition analysis method will be described in detail based on examples.
  • compositional analysis method of the present invention is not to be construed as being limited by the following specific examples.
  • Example 1 (Preparation of sample) Dioxane solution of three types of polymers: poly (methyl methacrylate), symbol “M”), polystyrene (symbol “S”), and polyethyl methacrylate (symbol “E”) (about 4% by mass, about 0.4 mg as polymer solid content) was cast on a sample pot and air-dried overnight to obtain a sample.
  • Table 4 shows the composition of each sample. Mass spectrometry was performed using LCMS-2020 manufactured by SHIMADZU as a detector, a DART ion source manufactured by IonSense, and an ion rocket manufactured by Biochromato as a sample heating unit.
  • FIG. 5 is a diagram showing the original spectrum (A) for "MSE_calib" and the reconstructed spectrum (B) after performing the above processing. It can be seen that both waveforms are similar. In other words, it was shown that the original spectral features were preserved by the above processing, making the data easier to use for compositional analysis.
  • FIG. 6 shows a part C ⁇ (MSE_Calib) of the C matrix after NMF decomposition and S T after the background and impurity spectra are removed by the CCA-filter.
  • Fig. 7 is the composition analysis result of the M, S, E ternary system.
  • hollow circles and "Calibrated" in the legend indicate that MSE_calib is used as a calibration sample, that is, the composition information of MSE_calib is given as known, and the intensity distribution matrix is intensity-corrected. This is the composition estimation result of
  • Non-calibrated is an environmental variable (i.e., 610 m / z: peak intensity of polydimethylsiloxane). It is the result of estimating the composition of a sample when the intensity distribution matrix is intensity-corrected.
  • FIG. 8 shows composition estimation results when the orthogonal constraint (1) is canceled from the NMF conditions.
  • the definition of each plot is the same as in FIG. From the results in FIG. 8, it can be seen that better estimation results are obtained when the NMF has orthogonality constraints.
  • FIG. 9 shows composition estimation results when only merging S is canceled from the MNF optimization conditions.
  • the definition of each plot is the same as in FIG. From the results of FIG. 9, it was found that better estimation results are obtained when the MNF optimization condition has merging S. Although the number of components can only be reduced to 36 by ARD, when the condition of Merging S is maintained, the composition can be estimated more accurately with the number of components of 9, which shows the importance of Merging S.
  • FIG. 10 shows composition estimation results (comparative example) when the CCA-filter is canceled.
  • the definition of each plot is the same as in FIG. In FIG. 10, the calculation result is that the sample "SE" contains about 1% of M even though the content of M is 0%. If 0% is estimated to be 1%, the estimation accuracy is not even several percent, which is not practical.
  • FIG. 11 shows the original spectrum of “MSE_calib” (upper, raw of MSE_calib) and the reconstructed spectrum (lower, reconstructed of MSE_calib) after the above processing with the CCA-filter removed. be. From this result, it was found that when the CCA-filter was canceled, neither DOX nor BG noise (see description of FIG. 5) could be removed.
  • Fig. 12 shows the composition estimation results when correction by environmental variables is not performed.
  • "Non-calibrated" represented by solid circles deviates greatly from Known, and it was found that the estimation results obtained are more accurate by correction using environmental variables.
  • Example 2 (Estimation based on samples not containing endmembers)
  • the samples described in Table 5 were prepared and the same test as in Example 1 was performed.
  • the symbols "M”, "S” and “E” mean the same as in Table 4 (Example 1).
  • each sample is adjusted so that the weight fraction of a single component (polymer) does not exceed 80%.
  • a two-dimensional mass spectrum was obtained for the specimen by the same method as in Example 1, and converted into a data matrix.
  • the first NMF processing already described, i.e. carried out. Note that estimation of the variance-covariance matrix of noise can be calculated using only X, so it is performed inside the module in the first NMF processing.
  • FIG. 15 illustrates C and BS as a result.
  • the composition estimated by this method is indicated by a circle, and the composition of the sample used (correct composition) is indicated by an asterisk. Both show good agreement.
  • the absence of plots at the vertices of C also indicates that the samples used did not contain endmembers, ie pure polymers of E, S and M.
  • BS in the figure shows the spectrum of the E, M, S pure polymer obtained by estimation and the spectrum of E, M, S actually measured, but there is no difference between the two, and the absolute intensity It was confirmed that it was well estimated, including
  • the pure polymer spectra of E, M, S and the composition of the mixture can be accurately estimated as shown in FIG. 15 BS. rice field. If there is test data, X test is continuously projected onto S and B to obtain C test .
  • Example 3 (Estimation from a specimen containing one endmember) Next, if one of the (training) specimens contains an endmember, the component content in the other specimen is arbitrary (it does not have to be the area outside the inscribed circle of the K ⁇ 1 dimensional simplex). We verified that.
  • Composition estimation was performed in the same manner as in Example 2, except that among the samples listed in Table 5, samples containing 40% by mass or more of S were excluded and M pure polymer was added as a sample.
  • FIG. 16 shows the results.
  • the solid circle points represent the estimation results of the sample used for estimating B and S, and the asterisks are the actual compositions. From this result, it was confirmed that S can be estimated with high accuracy, although there were no specimens outside the inscribed circle (hollow circles are specimens that were not used for estimation).
  • the non-orthogonality term is considered to be due to the incorporation of as a penalty term.
  • fragments of the material are directly pyrolyzed in the atmosphere without pretreatment, the resulting gas is analyzed by DART-MS, and unsupervised learning of the spectrum is performed to obtain a pure endmember spectrum of each component. was estimated, and a method was developed to represent the spectrum of the sample by linear summation of the sum. This allows the identification and quantification of the components of even materials of complex composition without the need for separation.
  • each peak intensity of the endmember spectrum can be estimated as an absolute intensity instead of a normalized relative value. It can be obtained as a ratio.
  • the present invention which breaks down the sample into components and finds the content ratio, can be a useful tool in material development, such as identifying key components that are factors in demonstrating the performance of materials. Furthermore, it is expected to have a wide range of applications, such as quality control such as content inspection of degraded components, legal compliance inspection such as hazardous material content, and inspection for detecting infringement of rights. Furthermore, it is expected that technologies such as deterioration and lifespan prediction of resins, detection of minute impurity components in resins and organic materials, quantitative analysis of smell and taste, traceability of food products, and tags to prevent falsification of origin can be expected to accelerate.
  • Plastics Europe proposed to the European Commission that plastic containers and packaging should contain 30% recycled plastic by 2030, promoting the use of recycled materials.
  • the present invention is expected to be applied to quality control, legal compliance inspection, etc., as a technique capable of accurately measuring the ratio of virgin agent and recycled material.

Landscapes

  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computing Systems (AREA)
  • Theoretical Computer Science (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Medical Treatment And Welfare Office Work (AREA)

Abstract

本発明の方法は、Ｋ種類から選択される成分を含む検体における、成分の含有量比を推定する方法であって、検体セットを加熱し、生ずるガス成分を順次イオン化し、マススペクトルを連続的に観測することと、マススペクトルから、検体ごとの二次元マススペクトルを得て、これらを複数まとめてデータ行列を得ることと、データ行列を非負値行列因子分解することと、基底スペクトル行列とデータ行列との正準相関分析によって強度分布行列を補正することと、補正後強度分布行列から特徴ベクトルを得て、ベクトル空間内で表現することと、Ｋ－１次元単体を設定し、エンドメンバーを決定することと、エンドメンバーと特徴ベクトルとから検体中の成分の含有量比を推定することと、を含み、アンビエントな条件で得られるマススペクトルからであっても、未知の混合物の成分の組成を推定することができる。

Description

検体に含まれる成分の含有量比の推定方法、組成推定装置、及び、プログラム
　本発明は、検体に含まれる成分の含有量比の推定方法、組成推定装置、及び、プログラムに関する。
　従来、多成分の混合系である検体における成分組成を推定するための方法としては、例えば、検体中の各成分を何らかの方法で分離して、それぞれの分画成分を同定、定量する方法が知られている。この同定方法の一つとして、質量分析装置によって得られるマススペクトルが使われることがあった。
　ところが、マススペクトルから得られる情報は、多数の相関関係を有する膨大な多次元データであるため、成分の同定につながる特徴的な信号を直感的、又は、経験的に抽出するのは困難であった。近年、非特許文献１に記載されているとおり、これらの多次元データを主成分分析によって次元削減し、特徴を探索する方法が用いられてきている。しかし、これら既知のデータ解析手法を用いても、ピーク強度がイオン化効率によって左右されるマススペクトルを用いて定量を行うことは本質的に困難であった。定量分析を行うためには、各化合物ごと有する未知変数であるイオン化効率をまず求める必要があり、そのためには検体量が既知、かつイオン化効率が同じとなる同位体検体を標準物質として用いる必要があった。
Ａｎａｌｙｔｉｃａｌ　ｃｈｅｍｉｓｔｒｙ，２０２０，ｖｏｌ．９２，ｉｓｓｕｅ２，Ｐ．１９２５－１９３３
　非特許文献１に記載された方法は、次元削減の過程で、もとの多次元データが有する多くの情報が失われ、成分の定性的な分析はできても、定量分析を行うことはできなかった。
　また、そもそも、成分の種類や試験環境等によって、成分ごとのイオン化効率が異なるという特性がある以上、得られるマススペクトルの大きさ（面積）の比は、混合物中における各成分の含有量の比を反映せず、これがマススペクトルによる定量分析の妨げの一つの要因となっていた。
　近年、質量分析方法の一つとして、アンビエントな条件（大気圧の開放空間にある）にある試料を前処理なしで質量分析するアンビエント脱離イオン化法（ＡＤＩ）が広く用いられるようになってきている。この方法は、未処理の試料や開放系の試料を固体のまま分析できるという利点がある一方で、各成分のイオン化率の不一致、及び、測定環境からの影響等によって、定量分析がより困難だと考えられてきた。
　上記事情を考慮し、本発明はアンビエントな条件で得られるマススペクトルからであっても、未知の混合物の成分の組成を推定することができる、成分の含有量比の推定方法を提供することを課題とする。
　また、本発明は、組成推定装置、及び、プログラムを提供することも課題とする。
　本発明者らは、上記課題を達成すべく鋭意検討した結果、以下の構成により上記課題を達成することができることを見出した。
　［１］　Ｋを１以上の整数としたとき、Ｋ種類の成分から選択される少なくとも１種の上記成分を含む推定対象検体における、上記成分の含有量比を推定する方法であって、上記Ｋ種類の成分から選択される少なくとも１種の成分を含み、互いに組成の異なるＫ個以上の学習用検体、及び、上記成分を含まないバックグラウンド検体を準備することと、上記推定対象検体、上記学習用検体、及び、バックグラウンド検体を含む検体セットに含まれるそれぞれの検体を加熱しながら、熱脱着、及び／又は、熱分解によって生ずるガス成分を順次イオン化し、マススペクトルを連続的に観測することと、加熱温度ごとに得られた上記マススペクトルを各行に格納し、上記検体ごとの二次元マススペクトルを得て、上記二次元マススペクトルの少なくとも２つ以上をまとめてデータ行列に変換することと、上記データ行列を非負値行列因子分解し、規格化された基底スペクトル行列とその強度分布行列の積に分解するＮＭＦ処理を行うことと、上記基底スペクトル行列と上記データ行列との正準相関分析によって強度分布行列に含まれるノイズ成分を抽出し、上記ノイズ成分の影響を減ずるよう上記強度分布行列を補正して、補正後強度分布行列を得ることと、上記補正後強度分布行列を上記検体のそれぞれに対応する小行列に分割し、上記小行列を特徴ベクトルとして、ベクトル空間内で上記検体のそれぞれを表現することと、上記特徴ベクトルの全てを内包するＫ－１次元単体を設定し、上記Ｋ－１次元単体におけるＫ個のエンドメンバーを決定することと、上記Ｋ個のエンドメンバーと、上記推定対象検体の特徴ベクトルとの、ユークリッド距離をそれぞれ計算し、上記ユークリッド距離の比により、上記推定対象検体中の上記成分の含有量比を推定することと、を含み、上記Ｋが３以上の場合、上記Ｋ－１次元単体に内接する超球体の外側の領域のそれぞれに、上記学習用検体の特徴ベクトルの少なくとも１つが位置する、又は、上記学習用検体が上記エンドメンバーの少なくとも１つを含む、方法。
　［２］　上記Ｋが２以上の整数であり、上記推定対象検体が、上記成分の混合物である、［１］に記載の方法。
　［３］　上記学習用検体が、上記エンドメンバーを含む、［１］又は［２］に記載の方法。
　［４］　上記エンドメンバーの決定が、上記学習用検体に付された判別ラベルに基づき行われる、［３］に記載の方法。
　［５］　上記エンドメンバーの決定が、上記学習用検体の上記特徴ベクトルの頂点成分分析によって行われる、［３］に記載の方法。
　［６］　上記エンドメンバーの決定が、上記Ｋ－１次元単体が最小体積となるように頂点を定めるアルゴリズムによって実施される、［１］に記載の方法。
　［７］　上記エンドメンバーの決定が、上記補正後強度分布行列を非負値行列因子分解し、上記検体中のＫ種類の上記成分の重量分率を表す行列と、Ｋ種類の個別の上記成分ごとのフラグメント存在強度を表す行列との積に分解する第２のＮＭＦ処理により行われる、［６］に記載の方法。
　［８］　上記学習用検体が、上記エンドメンバーを含まない、［６］又は［７］に記載の方法。
　［９］　上記補正後強度分布行列を得ることは、更に、上記強度分布行列の強度補正を行うことを含む、［１］～［８］のいずれかに記載の方法。
　［１０］　上記検体セットが、更に、キャリブレーション検体を含み、上記キャリブレーション検体は、上記Ｋ種類の成分を全て含み、かつ、その組成が既知であり、上記強度補正は、上記強度分布行列を規格化することを含む、［９］に記載の方法。
　［１１］　上記強度補正は、上記強度分布行列の少なくとも一部を、上記検体セットのうち対応する上記検体の質量と環境変数の積で割り付けることを含み、
　上記環境変数は、上記観測の際における上記成分のイオン化効率への影響を表す変数である、［９］に記載の方法。
　［１２］　上記環境変数が、上記観測の際に、雰囲気中に所定量含まれる分子量５０～１５００の化合物、又は、上記検体セットのそれぞれに所定量含まれる分子量５０～５００の有機低分子化合物のマススペクトルのピークの合計値である、［１１］に記載の方法。
　［１３］　Ｋを１以上の整数としたとき、Ｋ種類の成分から選択される少なくとも１種の上記成分を含む推定対象検体における、上記成分の含有量比を推定する組成推定装置であって、上記Ｋ種類の成分から選択される少なくとも１種の成分を含み、互いに組成の異なるＫ個以上の学習用検体、上記成分を含まないバックグラウンド検体、及び、上記推定対象検体を含む検体セットに含まれるそれぞれの検体を加熱しながら、熱脱着、及び／又は、熱分解によって生ずるガス成分を順次イオン化し、マススペクトルを連続的に観測する質量分析装置と、上記観測されたマススペクトルを処理する情報処理装置と、を備え、上記情報処理装置は、加熱温度ごとに得られた上記マススペクトルを各行に格納し、上記検体ごとの二次元マススペクトルを得て、上記二次元マススペクトルの少なくとも２つ以上をまとめてデータ行列に変換する、データ行列作成部と、上記データ行列を非負値行列因子分解し、規格化された基底スペクトル行列とその強度分布行列の積に分解するＮＭＦ処理を行う、ＮＭＦ処理部と、上記基底スペクトル行列と上記データ行列との正準相関分析によって強度分布行列に含まれるノイズ成分を抽出し、上記ノイズ成分の影響を減ずるよう上記強度分布行列を補正して、補正後強度分布行列を作成する、補正処理部と、上記補正後強度分布行列を上記検体セットのそれぞれに対応する小行列に分割し、上記小行列を特徴ベクトルとして、ベクトル空間内で上記検体のそれぞれを表現する、ベクトル処理部と、上記特徴ベクトルの全てを内包するＫ－１次元単体を設定し、上記Ｋ－１次元単体におけるＫ個のエンドメンバーを決定する、エンドメンバー決定部と、上記Ｋ個のエンドメンバーと、上記推定対象検体の特徴ベクトルとの、ユークリッド距離をそれぞれ計算し、上記ユークリッド距離の比により、上記推定対象検体中の上記成分の含有量比を推定する、含有量比計算部と、を含み、上記Ｋが３以上の場合、上記Ｋ－１次元単体に内接する超球体の外側の領域のそれぞれに、上記学習用検体の特徴ベクトルの少なくとも１つが位置する、又は、上記学習用検体が上記エンドメンバーの少なくとも１つを含む、組成推定装置。
　［１４］　上記エンドメンバー決定部は、上記学習用検体に付された判別ラベルに基づいて上記エンドメンバーを決定する、［１３］に記載の組成推定装置。
　［１５］　上記エンドメンバー決定部は、上記Ｋ－１次元単体が最小体積となるように頂点を定めるアルゴリズムによって上記エンドメンバーを決定する、［１３］に記載の組成推定装置。
　［１６］　上記エンドメンバー決定部は、上記補正後強度分布行列を非負値行列因子分解し、上記検体中のＫ種類の上記成分の重量分率を表す行列と、Ｋ種類の個別の上記成分ごとのフラグメント存在強度を表す行列との積に分解する第２のＮＭＦ処理により上記エンドメンバーを決定する、［１５］に記載の組成推定装置。
　［１７］　上記補正処理部は更に、上記強度分布行列の強度補正を行う、［１５］に記載の組成推定装置。
　［１８］　上記検体セットが、更に、キャリブレーション検体を含み、上記キャリブレーション検体が、上記Ｋ種類の成分を全て含み、かつ、その組成が既知であり、上記強度補正が、上記強度分布行列を規格化することを含む、［１７］に記載の組成推定装置。
　［１９］　上記強度補正は、上記強度分布行列の少なくとも一部を、上記検体セットのうち対応する上記検体の質量と環境変数の積で割り付けることを含み、
　上記環境変数は、上記観測の際における上記成分のイオン化効率への影響を表す変数である、［１７］に記載の組成推定装置。
　［２０］　上記質量分析装置は減圧ポンプを更に備え、上記環境変数が、上記減圧ポンプの稼働によって生ずる物質のマススペクトルのピークの合計値である、［１９］に記載の組成推定装置。
　［２１］　Ｋを１以上の整数としたとき、Ｋ種類の成分から選択される少なくとも１種の上記成分を含む推定対象検体における、上記成分の含有量比を推定する組成推定装置において用いられるプログラムであって、上記Ｋ種類の成分から選択される少なくとも１種の成分を含み、互いに組成の異なるＫ個以上の学習用検体、上記成分を含まないバックグラウンド検体、及び、上記推定対象検体を含む検体セットに含まれるそれぞれの検体を加熱しながら、熱脱着、及び／又は、熱分解によって生ずるガス成分を順次イオン化し、マススペクトルを連続的に観測する質量分析装置と、上記観測されたマススペクトルを処理する情報処理装置と、を備え、上記質量分析装置によって、加熱温度ごとに取得された上記マススペクトルを各行に格納し、上記検体ごとの二次元マススペクトルを得て、上記二次元マススペクトルの少なくとも２つ以上をまとめてデータ行列に変換する、データ行列作成機能と、上記データ行列を非負値行列因子分解し、規格化された基底スペクトル行列とその強度分布行列の積に分解するＮＭＦ処理を行う、ＮＭＦ処理機能と、上記基底スペクトル行列と上記データ行列との正準相関分析によって強度分布行列に含まれるノイズ成分を抽出し、上記ノイズ成分の影響を減ずるよう上記強度分布行列を補正して、補正後強度分布行列を作成する、補正処理機能と、上記補正後強度分布行列を上記検体のそれぞれに対応する小行列に分割し、上記小行列を特徴ベクトルとして、ベクトル空間内で上記検体のそれぞれを表現する、ベクトル処理機能と、上記特徴ベクトルの全てを内包するＫ－１次元単体を設定し、上記Ｋ－１次元単体におけるＫ個のエンドメンバーを決定する、エンドメンバー決定機能と、上記Ｋ個のエンドメンバーと、上記推定対象検体の特徴ベクトルとの、ユークリッド距離をそれぞれ計算し、上記ユークリッド距離の比により、上記推定対象検体中の上記成分の含有量比を推定する、含有量比計算機能と、を含み、上記Ｋが３以上の場合、上記Ｋ－１次元単体に内接する超球体の外側の領域のそれぞれに、上記学習用検体の特徴ベクトルの少なくとも１つが位置する、又は、上記学習用検体が上記エンドメンバーの少なくとも１つを含む、プログラム。
　本発明の方法（以下「本方法」ともいう。）は、検体セットに含まれるそれぞれの検体を加熱しながら、熱脱着、及び／又は、熱分解によって生ずるガス成分を順次イオン化し、マススペクトルを連続的に観測し、得られた加熱温度ごとのマススペクトルを各行に格納した検体ごとの二次元マススペクトルを得て、これらをまとめたデータ行列を非負値行列因子分解（ＮＭＦ）によって、規格化された基底スペクトル行列とその強度分布行列の積に分解し、正準相関分析によって強度分布行列に含まれるノイズ成分を抽出して補正し、各検体の特徴ベクトルを含む補正後強度分布行列を得て、その後の定量計算を行う点に特徴点の一つがある。
　この正準相関分析によって強度分布行列を補正する処理を行うことによって、（詳細は後述するが）より容易に、かつ、高精度で定量分析が可能になった。
　また、Ｋが２以上の整数であり、推定対象となる検体（推定対象検体）がＫ種類の成分の混合物である場合、本方法によって、推定対象検体中におけるＫ種類の成分の混合比（組成）を推定することができる。
　また、学習用検体がエンドメンバーを含む場合、より正確、かつ、より容易に検体の組成を推定することができる。
　また、学習用検体にエンドメンバーの判別するための判別ラベルが付されている場合、そのラベルの情報からエンドメンバーを決定できるため、より容易に検体の組成を推定することができる。
　また、判別ラベルが付されていない場合であって、学習用検体がエンドメンバーを含んでいる場合には、学習用検体の特徴ベクトルの頂点成分分析によってエンドメンバーを決定でき、検体の組成が推定できる。
　また、学習用検体がエンドメンバーを含んでいない場合であっても、Ｋ－１次元単体が最小体積となるように頂点を定めるアルゴリズムを用いることによってエンドメンバーを決定することで、検体の組成が推定できる。学習用検体にエンドメンバーを必須としない点で優れている。
　また、補正後強度分布行列が、強度分布行列の強度補正後のものであると、推定の正確性がより高まる。
　また、検体セットが更に、キャリブレーション検体を含む場合であって、強度補正が強度分布行列を規格化することを含む場合、より簡便な手順で、より高精度に検体の組成が推定できる。
　また、強度分布行列における各検体に対応する部分を検体の質量と環境変数の積で割り付けることによれば、観測の際における各成分のイオン化効率による影響が補正され、検体セットにキャリブレーション検体が含まれない場合であっても、検体の組成が高精度に推定できる。
　また、上記割り付けに使用される環境変数が観測の際の雰囲気中に所定量含まれる分子量５０～１５００の化合物、又は、検体セットのそれぞれの検体に所定量含まれる分子量５０～５００の有機低分子化合物のマススペクトルのピークの合計値である場合、検体セットにキャリブレーション検体が含まれない場合であっても、より簡便な操作で、より正確な組成推定結果が得られる。
本発明の実施形態に係る組成推定方法のフローチャートである。 Ｋが３のときのＫ－１次元単体（２次元単体、三角形）を表すイメージ図である。 本発明の実施形態に係る組成推定装置のハードウェア構成図である。 本発明の実施形態に係る組成推定装置の機能ブロック図である。 「ＭＳＥ＿ｃａｌｉｂ」についてのオリジナルなスペクトル（Ａ）と、処理を行った後、再構築したスペクトル（Ｂ）を表す図である。 ＮＭＦ分解後のＣ行列の一部Ｃδ（ＭＳＥ＿Ｃａｌｉｂ）と、ＣＣＡ－Ｆｉｌｔｅｒによってバックグラウンド及び不純物のスペクトルを除いた後のＳＴである。 Ｍ（ポリメチルメタクリレート）、Ｓ（ポリスチレン）、Ｅ（ポリエチルメタクリレート）三成分系の組成解析結果である。 ＮＭＦの条件から、直交制約を解除した場合の組成推定結果である。 ＭＮＦ最適化の条件から、ｍｅｒｇｉｎｇ　Ｓのみ解除した場合の組成推定結果である。 ＣＣＡ－ｆｉｌｔｅｒを解除した場合の結果（比較例）の組成推定結果である。 「ＭＳＥ＿ｃａｌｉｂ」についてのオリジナルなスペクトル（上、ｒａｗ　ｏｆ　ＭＳＥ＿ｃａｌｉｂ）と、ＣＣＡ－ｆｉｌｔｅｒを解除した上記で処理を行った後、再構築したスペクトルを表す図である。 環境変数による補正を行わなかった場合の組成推定結果である。 Ｋが３のときのＫ－１次元単体（２次元単体、三角形）を表すイメージ図である。 本発明の第２実施形態におけるエンドメンバーの決定方法の仕組みの説明図である。 実施例２の推定結果を表す図である。 実施例３の推定結果を表す図である。
　以下、本発明について詳細に説明する。
　以下に記載する構成要件の説明は、本発明の代表的な実施形態に基づいてなされることがあるが、本発明はそのような実施形態に制限されるものではない。
　なお、本明細書において、「～」を用いて表される数値範囲は、「～」の前後に記載される数値を下限値及び上限値として含む範囲を意味する。
［用語の定義］
　本明細書において使用される用語について説明する。なお、以下に説明のない用語については、当業者の間で普通に理解される意味で使用される。
　本明細書において、「成分」とは、各検体を構成する物質を意味し、推定対象検体、学習用検体、及び、キャリブレーション検体は、成分の組合せによって構成される混合物を意味する。
　「推定対象検体」とは、組成を推定すべき検体で、含まれている成分の数、及び、各成分の含有量はいずれも未知であってよい。
　「学習用検体」とは、後述するＫ－１次元単体のＫ個のエンドメンバー決定に必要な検体を意味し、これらの検体には、Ｋ種類の成分から選択される少なくとも１種の成分が含まれ、それぞれ組成が異なっている。含まれる成分数は特に制限されず１～Ｋ種類であってよい。なお、学習用検体には、推定対象検体と同一の組成のものが含まれていてもよい。すなわち、推定対象検体と学習用検体とは同一であってもよいが、学習用検体はそれぞれ異なる。
「エンドメンバー」とは、Ｋ－１次元単体の頂点に対応するベクトルを意味し、このエンドメンバーは、Ｋ種類のうちの１成分のみからなる検体（学習用検体、又は、仮想の検体）の特徴ベクトルに対応する。
「バックグラウンド検体」とは、推定対象検体、及び、学習用検体とは異なる検体であって、上記成分のいずれをも含まない検体を意味する。
「キャリブレーション検体」とは、推定対象検体、学習用検体、及び、バックグラウンド検体とは異なる検体であって、Ｋ種類の成分を全て含み、その組成が明らかな検体を意味する。
［本方法の第１実施形態］
　図１は、本発明の実施形態に係る組成推定方法（以下、「本方法」ともいう。）のフローチャートである。
　まず、推定対象検体、互いに組成の異なるＫ個以上の学習用検体、及び、バックグラウンド検体（以上をまとめて「検体セット」という。）が準備される（ステップＳ１０１）。なお、検体セットは、キャリブレーション検体を更に含んでもよい。
　推定対象検体は、Ｋ種類の成分のうち少なくとも１種の成分を含む。推定対象検体が１種類の成分を含む場合、本方法によってその成分の純度を推定することができる。この場合、推定対象検体には、上記「成分」と、それ以外の不純物とが含まれる。
　Ｋが２以上の場合、推定対象検体は、成分の混合物である。本方法におけるＫは１以上であれば、上限は特に制限されない。典型的には、Ｋは、１０００以下が好ましい。
　なお、学習用検体の個数は成分数Ｋに関係して、Ｋ個以上である必要がある。Ｋ個以上であれば、学習用検体の数としては特に制限されないが、より多い場合、後述するＫ－１次元単体の頂点の決定の精度がより高くなりやすく、結果として、より正確な組成推定結果が得られやすい。
　すなわち、成分数Ｋに対する学習用検体の数の比（学習用検体の数／成分数）は、１以上であり、２以上が好ましく、３以上がより好ましく、５以上が更に好ましく、１０以上が特に好ましい。一方、上限値は特に制限されないが、一般に、１０００以下が好ましい。
　学習用検体は、互いに組成が異なる。組成は、含まれる成分の種類と数、及び、各成分の含有量によって決まる。組成が同一とは、含まれる成分の種類、数、及び、それぞれの含有量が同一であることを意味し、いずれか１つが異なれば、それは「互いに組成が異なる」状態となる。
　なお、本方法において、学習用検体は、互いに組成が異なっていれば、それ以外の情報、すなわち、含まれる成分の種類、成分数、及び、各含有量は、推定対象検体の組成推定には必要ない。すなわち、学習用検体の組成は、未知であってよい。
　推定対象検体と学習用検体との相違点は、学習用検体はそれぞれ組成が異なることが必要である点である。なお、学習用検体は、推定対象検体と組成が同一のものを含んでいてもよい。
　学習用検体は、エンドメンバーを含んでいてもよい。すなわち、学習用検体は、Ｋ種類の成分から選択される１種類の成分からなるもの（単味のもの、ピュア）を含んでもよい。学習用検体がエンドメンバーを含む場合、より容易に、高精度な組成推定が可能である。
　このとき、エンドメンバーである学習用検体は、判別ラベルを含んでもよい。学習用検体が判別ラベルを含む場合、エンドメンバーの決定がより容易になる点で好ましい。
　次に、推定対象検体をそれぞれ加熱して、熱脱着、及び／又は、熱分解によって生ずるガスを順次イオン化し、マススペクトルを取得する（ステップＳ１０２）。
　マススペクトルを取得する方法としては特に制限されないが、アンビエントな条件にある試料を前処理なしで質量分析する方法が好ましい。このようなイオン化、及び、質量分析方法としては、ＤＡＲＴ（Ｄｉｒｅｃｔ　Ａｎａｌｙｓｉｓ　ｉｎ　Ｒｅａｌ　Ｔｉｍｅ）イオン源と呼ばれるイオン源と質量分析計とを組合わせたＤＡＲＴ－ＭＳと呼ばれる質量分析装置が知られている。
　質量分析計としては特に制限されず、精密質量分析が可能なものが好ましく、四重極型、及び、飛行時間（ＴＯＦ）型等のいずれでもよい。
　マススペクトルの取得について、具体的な条件は特に制限されないが、非限定的な一例を挙げると、検体セットに含まれるそれぞれの検体を順次５０℃／ｍｉｎの昇温速度で加熱し、５０～５５０℃の温度範囲で生じる熱分解ガスに対して５０ｓｈｏｔ／ｍｉｎの間隔でヘリウムイオンを噴射しガスをイオン化させ、横軸にｍ／ｚ、縦軸に温度を有する二次元ＭＳスペクトルを得るという手順が挙げられる。
　次に、加熱温度ごとに得られたマススペクトルを各行に格納し、検体ごとの二次元マススペクトルを得て、この二次元マススペクトルの少なくとも２つ以上をまとめてデータ行列に変換する（ステップＳ１０３）。
　データ行列の取得方法を具体的に説明する。
　ステップＳ１０２では、所定の昇温間隔ごとに連続的にマススペクトルが取得される。これらのマススペクトルは、そのままデータ行列の作成に使用されてもよいが、所定の昇温温度範囲ごとに平均して使用されてもよい。所定の昇温温度範囲ごとに平均してマススペクトルを１本にまとめることで、データ量を圧縮することができる。このような昇温温度範囲としては、例えば、１０～３０℃程度が挙げられる。
　また、各スペクトルにおいて、ピーク強度は規格化されていてもよい。規格化の方法としては、例えば、ピーク強度の二乗和が１になるように規格化する方法が挙げられる。
　このようにすることで、ある検体について、一回の測定することで、加熱温度ごと（又は、所定の範囲ごとにまとめられた加熱温度ごと）に、所定の本数（まとめ方によって異なる、例えば、２０本等）のマススペクトルが得られる。
　このマススペクトルが各行に格納され、加熱温度が各列に格納されると、検体ごとの二次元マススペクトルが得られる。
　こうして、検体ごとの二次元マススペクトルが得られたら、これらの少なくとも２つ以上をまとめてデータ行列Ｘに変換する。
　データ行列Ｘの作成に用いられる二次元マススペクトルの個数は２つ以上であれば特に制限されないが、検体セットに含まれる全検体の二次元マススペクトルが用いられることが好ましい。
　なお、１つの検体について２回以上、測定が行われる場合には、２回以上の測定で得られた二次元マススペクトルの一部、又は、全てがデータ行列Ｘの作成に用いられてもよい。
　次に、得られたデータ行列を非負値行列因子分解（ＮＭＦ）し、規格化された基底スペクトル行列と、その強度分布行列との積に分解する（ステップＳ１０４）。
　非負値行列因子分解について説明する。
　スペクトルの強度とサンプルの成分組成の間に線形加成性が成り立つという仮定の下で、スペクトルのデータ行列は一般に、
Figure JPOXMLDOC01-appb-M000001
と行列分解できる。ここでデータ行列Ｘはチャンネル数ＭのスペクトルをＮ本積み重ねてできるＮ×Ｍのデータ行列、Ｓは基底スペクトル行列、Ｃはその強度分布行列である。データセットに含まれるすべての成分数をＫ１とすればＣはＮ×Ｋ１行列、ＳはＭ×Ｋ１行列になる。実際の計測にはノイズが含まれるので、
Figure JPOXMLDOC01-appb-M000002
と誤差行列Ｎを考慮にいれ、Ｎを最小化するようにＣ，Ｓ行列を定める。行列右肩のＴは転置行列を意味する。ここでガウスノイズを仮定すれば
Figure JPOXMLDOC01-appb-M000003
なるＣ、Ｓについての最適化問題として定式される。なお、
Figure JPOXMLDOC01-appb-M000004
は、行列のフロベニウスノルム、
Figure JPOXMLDOC01-appb-M000005
は、ベクトルのＬ１ノルム、
Figure JPOXMLDOC01-appb-M000006
は、ベクトルのＬ２ノルムを表すものとする。
　なお、上式におけるＣ、Ｓは唯一解が定まらないので、物理的要請からスペクトルが満たすであろう拘束条件を課す。ＮＭＦは、組成比やスペクトル強度が非負の値をとるという自明な要請から、
Figure JPOXMLDOC01-appb-M000007
  
として立式される。しかし、非負条件だけでは唯一解が定まらないので、測定装置に応じたさらなる拘束条件を適用することが求められる。以下では、本ステップにおけるＮＭＦの拘束条件について説明する。
（ＮＭＦ更新式について）
　本ステップにおけるＮＭＦ処理は、以下に述べる制約を組み込みやすいＨｉｅｒａｒｃｈｉｃａｌ　ａｌｔｅｒｎａｔｉｎｇ　ｌｅａｓｔ　ｓｑｕａｒｅ　（ＨＡＬＳ）により、Ｃ、Ｓの列ベクトルワイズに交互に更新していくことで、上記式（数７）をＣ、Ｓの各列についての凸最適化問題として解くことが好ましい。
　本方法では、各検体を加熱しながら、熱脱着、及び／又は、熱分解させて生じるガスについて質量分析を行うため、熱分解後の成分数Ｋ１を事前情報として得ることは難しく、データから推定する必要がある。推定方法としては、Ａｕｔｏｍａｔｉｃ　ｒｅｌｅｖａｎｃｅ　ｄｅｔｅｒｍｉｎａｔｉｏｎ　（ＡＲＤ）が好ましい。これは強度分布行列Ｃにスパース性を仮定するもので、ある成分の熱分解はその成分を含む検体の特定温度でのみ起こることを考えれば、自然な仮定といえる。これによりなるべく少ない成分数にする力と上記式（数７）を最小化する力が釣り合い、最適な成分数を求めることができる。
（直交制約）
　質量分析計は、一般に０．１～０．００１ｍ／ｚ程度の分解能を有する。したがって、異なるガス成分のｍ／ｚが小数点以下まで一致し、同一のチャネルを占めるという状況は稀にしか生じない。したがって、基底スペクトル行列Ｓの列ベクトル間にはある程度の直交性を仮定するのが妥当である。
　これには良い局所解に導く効果のみならず、強度の強いピークにのみフィッティングが集中することを防ぐ効果も期待できる。一方で、完全な直交かつ非負は強すぎる制約となってしまうので、ラグランジュ未定乗数を有するペナルティー項を加える形でソフトに直交制約を課すことが好ましい（ソフト直交制約）。
（成分のマージング）
　ＣとＳの更新途中で、ほとんど同じ基底スペクトルが複数得られた場合は一つのスペクトルにまとめることで正しい成分数ｎを得ることができる。ＥＡＬＳのスペクトライメージをＮＭＦに適用すると、Ｓのｋ番目の列ベクトルをＳ：ｋと書いて、ＳＴ ：ｋＳ：ｍ＞ｔｈｒｅｓｈｏｌｄなるｋ＜ｍが存在するときに、
Figure JPOXMLDOC01-appb-M000008
によって一本化できる（Ｍｅｒｇｉｎｇ　Ｃ）。ここでＴｈｒｅｓｈｏｌｄは０．９～０．９９で与えられるハイパーパラメータである。これと同様に類似Ｓのマージングを行い、更に、Ｃの複数の列ベクトルが平行であったときも、対応するＳの列ベクトルをＣの重みに応じて重ね合わせることが好ましい。すなわち、
Figure JPOXMLDOC01-appb-M000009
なるｋ≠ｍが存在するときに、
Figure JPOXMLDOC01-appb-M000010
と一本化する（Ｍｅｒｇｉｎｇ　Ｓ）。これには同位体ピークやオリゴマーピークを一つの基底スペクトルにまとめる効果があると考えられる。
（更新式の立式）
　上記によって得られる更新式について説明する。
　モデルの全同時分布は
Figure JPOXMLDOC01-appb-M000011
　ここではσ２は分散パラメータ、
Figure JPOXMLDOC01-appb-M000012
の要素λｋは成分ｋがデータ中に占める重要度を表す逆ガンマ分布に従うパラメータである。ノイズにガウス分布を仮定したので
Figure JPOXMLDOC01-appb-M000013
　Ｃが従う指数分布は、Ｋを初期値として入力する十分大きな任意の成分数（５０程度）として、
Figure JPOXMLDOC01-appb-M000014
  
　λｋが従う逆ガンマ分布は、
Figure JPOXMLDOC01-appb-M000015
　ここでａはａ－１が極めて小さな正の数となるようなハイパーパラメータ、ｂはデータ行列の平均値μＸにより
Figure JPOXMLDOC01-appb-M000016
  
として定まるパラメータである。また、Ｃの要素の期待値は
Figure JPOXMLDOC01-appb-M000017
ｐ（Ｘ，Ｃ，Ｓ）の負の対数尤度は
Figure JPOXMLDOC01-appb-M000018
　ここで、ガウス分布を仮定したときの分散パラメータσ2は、
Figure JPOXMLDOC01-appb-M000019
で与えられ、また、λについて凸関数なので
Figure JPOXMLDOC01-appb-M000020
を得ることができる。
　次に、Ｓの更新式について説明する。まず、負対数尤度関数からＳ：ｋに関連する項Ｊ（Ｓ：ｋ）のみを取り出すことを考える。Ｘに対するｋ成分の寄与、
Figure JPOXMLDOC01-appb-M000021
を用いて、
Figure JPOXMLDOC01-appb-M000022
と書けるので、定数項を除いて、
Figure JPOXMLDOC01-appb-M000023
　上記は、Ｓｈｉｇａ，Ｍ．ｅｔ　ａｌ．Ｓｐａｒｓｅ　ｍｏｄｅｌｉｎｇ　ｏｆ　ＥＥＬＳ　ａｎｄ　ＥＤＸ　ｓｐｅｃｔｒａｌ　ｉｍａｇｉｎｇ　ｄａｔａ　ｂｙ　ｎｏｎｎｅｇａｔｉｖｅ　ｍａｔｒｉｘ　ｆａｃｔｏｒｉｚａｔｉｏｎ．　Ｕｌｔｒａｍｉｃｒｏｓｃｏｐｙ　１７０，４３－５９（２０１６）．をベースに導出された。ここで、上述したＳに加える直交制約によるペナルティー項を加える。
Figure JPOXMLDOC01-appb-M000024
　ここで、ξｋはラグランジュ未定乗数、
Figure JPOXMLDOC01-appb-M000025
  
である。ｓ：ｋで微分して
Figure JPOXMLDOC01-appb-M000026
Figure JPOXMLDOC01-appb-M000027
とおいて、ｓ：ｋについての更新式
Figure JPOXMLDOC01-appb-M000028
を得る。完全直交条件を与えるξｋは、
Figure JPOXMLDOC01-appb-M000029
より、
Figure JPOXMLDOC01-appb-M000030
　より見積もることができる。ソフト直交条件はξｋに［０，１］のハイパーパラメータｗをかけて、完全直交と直交制約なしとの間で自由に直交度合いを調整することができる。最終的な更新式は、
Figure JPOXMLDOC01-appb-M000031
　次に、Ｃの更新式について説明する。負対数尤度関数からＣ：ｋに関連する項Ｊ（Ｃ：ｋ）のみを取り出すことを考える。
Figure JPOXMLDOC01-appb-M000032
Figure JPOXMLDOC01-appb-M000033
と規格化しておけば、更新式は
Figure JPOXMLDOC01-appb-M000034
より、
Figure JPOXMLDOC01-appb-M000035
　具体的なアルゴリズムは以下のとおりである。
Figure JPOXMLDOC01-appb-M000036
　図１に戻り、次に、基底スペクトル行列とデータ行列とを正準相関分析して補正後強度分布行列を得る（ステップＳ１０５）。
　具体的には、基底スペクトル行列とデータ行列との正準相関分析によって強度分布行列に含まれるノイズ成分を抽出し、ノイズ成分の影響を減ずるよう強度分布行列を補正して、補正後強度分布行列を得る。
　ＮＭＦはデータ行列の低ランク近似であるため、実際にはｉ番目のスペクトルには存在しない成分ｋであっても、存在するとしたほうが最小二乗の意味で近似が良くなるのであればＣｉｋ＞０とする。多くの場合このようなＣｉｋは非常に小さく、ＮＭＦ解析において問題にならない場合が多い。
　しかし、微量成分検出においては、ｊ番目のスペクトルに微量実在するＣｊｋ＞０と、ＮＭＦアーティファクトであるＣｉｋ＞０を区別し、ゴーストピークであるＣｉｋには０を代入することが好ましい。ＮＭＦアルゴリズムのアーティファクトに由来する偽ピークをノイズの１つとして除くことによりより精度の高い推定結果が得られるからである。
　本ステップでは、上記を解決するために、正準相関分析（Ｃａｎｏｎｉｃａｌ　Ｃｏｒｒｅｌａｔｉｏｎ　Ａｎａｌｙｓｉｓ　）を用いる。本ステップは本方法の特徴点の一つであり、本発明者らは本ステップを、Ｃａｎｏｎｉｃａｌ　Ｃｏｒｒｅｌａｔｉｏｎ　Ａｎａｌｙｓｉｓ　（ＣＣＡ）　Ｆｉｌｔｅｒと命名した。
　ＣＣＡ　ｆｉｌｔｅｒは、概念的には、ＮＭＦから出力される基底スペクトルＳＴの各成分が実際に元データに含まれていたかをサンプルワイズにスキャンし、元データに類似のピークパターンが見られなければそのサンプルのスペクトルから消去することによる。Ｌ番目のサンプルから生成されるスペクトルがデータ行列Ｘのｌ行からｌ′行までに格納されていた場合、その部分を抜き出した小行列をＸδ（Ｌ)と書く。すなわちδ（Ｌ）＝｛ｌ，ｌ＋１，…，ｌ′－２，ｌ′－１｝ある。
　ＣＣＡは従来、二系統のデータ間の相関を調べる手法として知られており、以下で定式される。
　N個のｐ次元ベクトルからなる行列を
Figure JPOXMLDOC01-appb-M000037
Ｎ個のｑ次元ベクトルデータからなる行列を
Figure JPOXMLDOC01-appb-M000038
として、
Figure JPOXMLDOC01-appb-M000039
ここで、
Figure JPOXMLDOC01-appb-M000040
である。
　すなわち、ある軸にＹとＺを投影して、その相関係数ρが最大となるような解（ａＴ，ｂＴ）Ｔを見つけることがＣＣＡの問題設定である。
　（ａＴ，ｂＴ）Ｔと相関係数ρは次式の一般固有値問題の固有ベクトルと固有値として与えられ、すべての固有ベクトル（ａＴ，ｂＴ）Ｔをまとめて（ＡＴ，ＢＴ）Ｔと書くと、
Figure JPOXMLDOC01-appb-M000041
　いま、ある成分ｋの基底スペクトルＳ：ｋが、サンプルＬに含まれているか、いないかを判定したいとする。まず、データ行列Ｘから該当の部分データ行列Ｘδ（Ｌ）を抜き出す。一方ＳＴは、Ｓ：ｋ Ｔと似たすべてのスペクトル群をＳｓｉｍ Ｔ、それとは大きく異なるスペクトル群をＳｏｕｔ Ｔと、二つの行列に分け、
Figure JPOXMLDOC01-appb-M000042
としてＣＣＡを行う。ここである値ｔ∈（０，１］に対して、ρ＞ｔ１なるすべての固有値を取り出し、それに対応する固有ベクトルｂの第一成分がその最大成分に対してｂ１＞ｔ２ｂｍａｘを満たせば、成分ｋの基底スペクトルＳ：ｋがサンプルＬに含まれていると判定する。そうでなければ、対応するＣδ（Ｌ）ｋ＝０とする。
　ここで、ｔ１、ｔ２はハイパーパラメータで、ＮＭＦのフィッティング精度、すなわちＮＭＦハイパーパラメータであるｔｈｒｅｓｈｏｌｄの影響を強く受けるので、慎重に吟味する必要がある。経験的にｔ１＝０．９ｔｈｒｅｓｈｏｌｄ，ｔ２∈［０．１，０．３］程度がよい結果を与える。
　上記は、ＣＣＡ　ｆｉｌｔｅｒをＮＭＦアーティファクト（ノイズ）に適用する例について説明したが、上記のノイズ除去は、ＮＭＦアーティファクトだけでなく、バックグラウンドやコンタミネーションの除去にも適用できる。
　すなわち、ＣＣＡ　ｆｉｌｔｅｒで除去できるノイズとしては、例えば、（１）ＮＭＦのアーティファクト（いわゆる偽ピーク）、（２）装置の周辺環境から紛れ込む微量化学物質由来のピーク（バックグラウンド）、（３）検体中に含まれるが、分析に含めるとその本質をゆがめる恐れのある物質由来のピーク（コンタミネーション）が挙げられる。
　本方法で扱う検体セットには、バックグランド検体が含まれているため、
Figure JPOXMLDOC01-appb-M000043
　とすれば、基底スペクトルＳ：ｋがバックグランドファイルに含まれていたどうかが判定できる。もし含まれていたのであれば、それは環境由来のシステムノイズであるから、Ｃ：ｋ＝０とすればサンプルから当該スペクトルを取り除くことができる。
　なお、更に、バックグラウンド検体がコンタミネーション物質を含む場合、上記と同様の手順によって、コンタミネーション由来のノイズも除去することができる。コンタミネーション物質としては特に制限されないが、例えば、検体の作製やロードに用いた溶媒等が挙げられる。
　本ステップは更に、強度分布行列の強度補正を行うステップを含んでいてもよい。
　例えば、ＤＡＲＴ－ＭＳは、空気中の水やアンモニアをヘリウムイオンで一次的にイオン化し、これによりサンプルを二次的にイオン化する。そのため、湿度等、大気条件によってピーク強度が変化する場合がある。またロードする検体質量も検体ごとに異なっているのが普通であり、これもピーク強度に影響を及ぼす。また、検体に含まれる各成分のイオン化効率も一般にはそれぞれ異なり、これもピーク強度に影響を及ぼす。
　強度分布行列の強度補正によれば、上記の影響を減じて、より正確な定量結果を得ることができる。
　全検体が、Ｋ種類の成分を全て含み、かつ、その組成が既知であるキャリブレーション検体を含む場合、上記強度補正は、上記強度分布行列を規格化することによって実施できる。
　強度分布行列を規格化するとは、Ｃ行列を検体ごとに規格化することであり、より具体的には、検体ごとにすべての要素を足して１になるようにする操作を意味する。
　Ｃ行列をサンプルワイズに分割し、サンプルＬに由来するＣ行列をＣδ（Ｌ）とかくと、
Figure JPOXMLDOC01-appb-M000044
組成分析においては分解温度の情報は不要であるので、
Figure JPOXMLDOC01-appb-M000045
を縦方向に加算して得られるベクトル
Figure JPOXMLDOC01-appb-M000046
をサンプルＬの特徴量ベクトルとする。次にエンドメンバーとして測定したｐ個の単成分サンプルの特徴ベクトルを並べて、エンドメンバー行列
Figure JPOXMLDOC01-appb-M000047
を作る。混合サンプルの特徴ベクトル
Figure JPOXMLDOC01-appb-M000048
はエンドメンバーの特徴ベクトルの線形和で書けるので、その係数ベクトル
Figure JPOXMLDOC01-appb-M000049
を使って、
Figure JPOXMLDOC01-appb-M000050
と書ける。
　一つのキャリブレーション検体は必ず一つ以上の熱分解サンプルを与えるので、熱分解後の成分数はＰは列フルランクな縦長の行列である。したがって一般的には解は不能で、最小二乗の意味での近似解は二次計画法で解くことができる。
　本方法において得られるマススペクトルでは、イオン化効率が化合物ごとに異なるうえ、熱分解効率も異なり、化合物量とピーク強度は比例しない。
　しかし、ＮＭＦによって成分ｋのイオン化効率を含めたピーク強度比はＳ：ｋ Ｔの中にすでに反映されており、スペクトル間の強度比はＣ：ｋに反映されている。
　このＮＭＦによって形成された階層的データ構造のおかげで、エンドメンバーにおけるピーク強度比さえ考慮にいれれば、未知のイオン化効率を求めなくても、より正確な定量的な組成分析が可能になる。
　一方、検体セットがキャリブレーション検体を含まない場合、強度分布行列の少なくとも一部（好ましくは全部）を、対応する検体の質量と環境変数の積で割り付けることによって、より正確な定量分析が可能となる。
　ここで、環境変数とは、観測の際における各成分のイオン化効率への影響を表す変数である。
　イオン化効率が高い成分であれば、ピーク強度は大きくなり、イオン化効率が低い成分であれば、ピーク強度が小さくなる。
　また、検体の質量についても、これが大きくなれば、一般的にピーク強度は大きく成り、小さく成れば、ピーク強度は小さくなる。
　検体セットがキャリブレーション検体を含まない場合、強度分布行列を規格化せず、強度分布行列が保持している検体間の強度情報を用いることができるよう、Ｃ行列を次のように補正する。
Figure JPOXMLDOC01-appb-M000051
ここでｍＬはサンプルＬの質量、ＩＬは６１０ｍ／ｚに全温度域で観測される、ＤＡＲＴ－ＭＳ付属のオイルポンプから生じる（ポリ）ジメチルシロキサンのピーク強度である。ＩＬはその測定時の待機状態を反映しており、実際雨の日と晴れの日ではピーク強度に３倍程度の差がある。これにより補正した
Figure JPOXMLDOC01-appb-M000052
縦方向に加算して得られるベクトル
Figure JPOXMLDOC01-appb-M000053
をサンプルＬの特徴ベクトルとする。あとは同様にエンドメンバーの特徴ベクトルからエンドメンバー行列
Figure JPOXMLDOC01-appb-M000054
をつくり、
Figure JPOXMLDOC01-appb-M000055
よりｒＬを求める。
　なお、上記ではＩＬは、（ポリ）ジメチルシロキサン由来の６１０ｍ／ｚとしたが、上記に制限されず、観測の際に雰囲気中に所定量に含まれる分子量５０～１５００の化合物に由来するもの、及び／又は、全検体のそれぞれに所定量含まれる分子量５０～５００の有機低分子化合物に由来するものであることが好ましい。
　図１に戻り、次に、補正後強度分布行列を検体セットに含まれる検体のそれぞれの検体に対応する小行列に分割し、これを特徴ベクトルとしてベクトル空間内で各検体を表現する（ステップＳ１０６）。
　さらに、これらの特徴ベクトル（それぞれ各検体に対応する）の全てを内包するＫ－１次元単体を設定し、Ｋ－１次元単体におけるＫ個のエンドメンバーを決定する（ステップＳ１０７）。
　学習用検体がエンドメンバーを含む場合であって、学習用検体に判別ラベル（当該学習用検体が「エンドメンバー」であると判別するためのラベル）が付されている場合、そのラベルを元にエンドメンバーが決定される。
　学習用検体がエンドメンバーを含む場合であって、学習用検体に判別ラベルが付されていない場合、すなわち、学習用検体にエンドメンバーが含まれてはいるが、どの学習用検体がエンドメンバーかが不明である場合、特徴ベクトルの頂点成分分析（Ｖｅｒｔｅｘ　Ｃｏｍｐｏｎｅｎｔ　Ａｎａｌｙｓｉｓ）によってエンドメンバーが決定される。
　学習用検体がエンドメンバーを含まない場合には、Ｋ－１次元単体が最小体積となるように頂点を定めるアルゴリズムによってエンドメンバーが決定される。なお、Ｋ－１次元単体が最小体積になる、とは、２次元単体（３角形）の場合は、最小面積を意味する。
　本方法の一形態においては、学習用検体がエンドメンバーを含まない場合でも、上記Ｋ－１次元単体に内接する超球体の外側の領域のそれぞれに、学習用検体の特徴ベクトルの少なくとも１つが位置しているため（Ｋが３以上の場合）、精度よくエンドメンバーを決定することができる。
　図２は、Ｋが３のときのＫ－１次元単体（２次元単体、三角形）を表すイメージ図である。Ｋ－１次元単体１０は、学習用検体１６、１７、及び、１８によって決定されたエンドメンバー１３、１４、及び、１５を頂点とする三角形である。
　図２のように、学習用検体１６、１７、及び、１８が、Ｋ－１次元単体の内接超球１２（この場合「円」）の外側の領域１９（ハッチング）のそれぞれに位置していると、精度の高い定量分析結果が得られる。
　なお、上記はＲｏｂｕｓｔ　Ｖｏｌｕｍｅ　Ｍｉｎｉｍｉｚａｔｉｏｎ－Ｂａｓｅｄ　Ｍａｔｒｉｘ　Ｆａｃｔｏｒｉｚａｔｉｏｎ　ｆｏｒ　Ｒｅｍｏｔｅ　Ｓｅｎｓｉｎｇ　ａｎｄ　Ｄｏｃｕｍｅｎｔ　Ｃｌｕｓｔｅｒｉｎｇ，　ＩＥＥＥ　ＴＲＡＮＳＡＣＴＩＯＮＳ　ＯＮ　ＳＩＧＮＡＬ　ＰＲＯＣＥＳＳＩＮＧ，ＶＯＬ．　６４，ＮＯ．　２３，ＤＥＣＥＭＢＥＲ　１，２０１６をベースに考案された。
　Ｋ－１次元単体が最小体積となるように頂点を定めるアルゴリズムとしては特に制限されず、公知のアルゴリズムが適用できる。このようなアルゴリズムとしては、例えば、Ｍｉｎｉｍｕｍ　ｖｏｌｕｍｅ　ｃｏｎｓｔｒａｉｎｔ－ＮＭＦ法、Ｍｉｎｉｍｕｍ　Ｄｉｓｔａｎｃｅ　Ｃｏｎｓｔｒａｉｎｔ－ＮＭＦ法、Ｓｉｍｐｌｅｘ　ｉｄｅｎｔｉｆｉｃａｔｉｏｎ　ｖｉａ　ｓｐｌｉｔ　ａｕｇｍｅｎｔｅｄ　Ｌａｇｒａｎｇｉａｎ法、及び、Ｓｉｍｐｌｅｘ　ｖｏｌｕｍｅ　ｍｉｎｉｍｉｚａｔｉｏｎ法等が挙げられる。
　次に、Ｋ個のエンドメンバーと、推定対象検体の特徴ベクトルとのユークリッド距離をそれぞれ計算し、検体中の各成分の含有量比を推定する（ステップＳ１０８）。
　なお、本方法の他の一形態として、学習用検体が、エンドメンバーを少なくとも１つ含む（学習用検体の少なくとも１つがエンドメンバーである）場合、他の学習用検体の特徴ベクトルは、Ｋ－１次元単体に内接する超球体の内側の領域に位置していてもよい。図１３は、図２と同様に、Ｋが３のときのＫ－１次元単体（２次元単体、三角形）を表すイメージ図である。Ｋ－１次元単体６０は、エンドメンバーである学習用検体６１、それ以外の学習用検体６２、及び、学習用検体６３によって決定されたエンドメンバー６４、及び、６５を頂点とする三角形である。図２と異なるのは、学習用検体６２、及び、６３が、Ｋ－１次元単体に内接する超球体の内側の領域に位置していることである。
　なお、学習用検体の少なくとも１つがエンドメンバーである場合であっても、他の学習用検体は、Ｋ－１次元単体に内接する超球体上に位置していてもよいし、外側の領域に位置していてもよい。すなわち、この場合、他の学習用検体の位置は任意でよい。
　なかでも、より優れた本発明の効果が得られる点では、他の学習用検体は、エンドメンバーである学習用検体とは異なるエンドメンバーの成分（他のエンドメンバーの成分）を２０質量％以上含むことが好ましく、４０質量％以上含むことがより好ましい。
　本方法によれば、溶媒に溶解しにくく、従来の組成分析が行いにくい固体の検体であっても、特段の前処理なしに分析に供することができ、かつ、正確な定量分析が可能になる。本方法は、特に、複数の成分の混合物である検体の品質管理、及び、不良原因の調査等に適用した場合、結果を得るために必要な時間を大幅に短縮することができる。
［本方法の第２実施形態］
　本方法の第２実施形態は、すでに説明した第１実施形態におけるステップＳ１０７において、学習用検体がエンドメンバーを含まない場合、すなわち、Ｋ－１次元単体頂点に位置する学習用検体が１つもない場合のエンドメンバー（頂点）の決定方法を含む。
　本実施形態に係る方法の特徴点の一つは、エンドメンバーの決定が、検体ごとのフラグメント存在強度を表す行列を非負値行列因子分解し、検体中のＫ種類の成分の重量分率を表す行列と、Ｋ種類の成分の各単体ごとのフラグメント存在強度を表す行列との積に分解する第２のＮＭＦ処理により行われる点にある。
　図１４は、その仕組みを説明するための模式図である。
　各検体は、Ｋ種類の成分の混合物（以下では、「ポリマー」の混合物として説明する）であるが、実際に観測できるのは、各ポリマーが熱分解して生じたＭ種類のフラグメントガスのスペクトルSの適当な線形和である。
　本実施形態に係る上記特徴点は、その係数行列を
Figure JPOXMLDOC01-appb-M000056
 とすれば、
Figure JPOXMLDOC01-appb-M000057
 であり、仮にＫ個のポリマーの単体について測定できたとして、それが
Figure JPOXMLDOC01-appb-M000058
 と観測されたとすると、
Figure JPOXMLDOC01-appb-M000059
 と、二度目のＮＭＦ処理を行えば、組成Ｃが得られる、という着想に基づいている。
　なお、ここで
Figure JPOXMLDOC01-appb-M000060
 はＫ個の単体（ポリマー）から生じるＭ個のフラグメントの存在量を表す。
　以下では、図１のフローにおける、ステップＳ１０４のＮＭＦ処理（以下、「１回目のＮＭＦ処理」ともいう。）と、ステップＳ１０７におけるエンドメンバーの決定において実施されるＮＭＦ処理（以下「２回目のＮＭＦ処理」ともいう。）の違い等について詳述する。なお、以下に説明のないステップ等については、原則として、第１実施形態における各ステップと同様である。
　以下の説明において、第１実施形態における説明で使用した各記号等を以下の表のとおり置き換えて説明する。
Figure JPOXMLDOC01-appb-T000061
　また、以下の説明において、多くは信号処理で一般的に使われている表記法に従う。
Figure JPOXMLDOC01-appb-M000062
は、は非負、又は、実数のＮ×Ｍ行列を表す。
　ある行列
Figure JPOXMLDOC01-appb-M000063
 に対して、
Figure JPOXMLDOC01-appb-M000064
 はそれぞれn番目の行ベクトル、m番目の列ベクトル、(n, m)-要素を表す。ＸＴはＸの転置行列、
Figure JPOXMLDOC01-appb-M000065
 はフロベニウスノルムを表す。
Figure JPOXMLDOC01-appb-M000066
 はＸのｎ行目の
Figure JPOXMLDOC01-appb-M000067
 を表す。Ｘが正方行列のとき、Ｔｒ（Ｘ）はトレース、ｄｉａｇ（Ｘ）は対角成分からなる対角行列を表す。１Ｎと１１Ｎはすべての要素が１のＮ次元ベクトル又は（Ｎ，Ｎ）行列を表す。ＩＮはＮ次元の単位行列を表す。
（１回目のＮＭＦ）
Figure JPOXMLDOC01-appb-M000068
　入力:データ行列
Figure JPOXMLDOC01-appb-M000069
出力：分布行列
Figure JPOXMLDOC01-appb-M000070
 基底フラグメントスペクトル
Figure JPOXMLDOC01-appb-M000071
 　ここで、Ｎ：サンプル数、Ｔ：温度帯域分割数、Ｄ：チャネル数、Ｍ：基底スペクトル数である。
　一回目のＮＭＦは、志賀らによって提案されたＡＲＤ－ＳＯ－ＮＭＦに対して以下３点を主な変更点として加え、ＭＳ（Ｍａｓｓ　ｓｐｅｃｔｒｏｍｅｔｒｙ）のデータ解釈によく適合するように開発した。
１）チャネルごとのガウスノイズの分散共分散行列
Figure JPOXMLDOC01-appb-M000072
 を自然同位体ピークに基づき見積もること
２）ソフト直交制約を基底フラグメントスペクトル間に適用すること
３）似た強度分布を有するフラグメントスペクトル同士を統合すること（マージング条件の拡張）
　１）については、以下しばらくの間、分散σ２で独立等分散（ｉ．ｉ．ｄ）のガウスノイズを仮定し、途中で分散共分散行列Ｒを導入する。ノイズにｉ．ｉ．ｄを仮定すると、データ行列Ｘの確率生成モデルは、
Figure JPOXMLDOC01-appb-M000073
 と書ける。基底スペクトル数Ｍは未知なので、自動的に推定できるように、ａｕｔｏｍａｔｉｃ　ｒｅｌｅｖａｎｃｅ　ｄｅｔｅｒｍｉｎａｔｉｏｎ　（ＡＲＤ）を分布行列Ａのスパース性に基づいて導入する。まず、基底成分ごとにλｍでパラメトライズされた指数分布をＡの事前分布として仮定する。
　すなわち、
Figure JPOXMLDOC01-appb-M000074
 に対して、
Figure JPOXMLDOC01-appb-M000075
  
全体の確率モデルとしては、
Figure JPOXMLDOC01-appb-M000076
 と書ける。ここで、ｐ（Ｓ）は、
Figure JPOXMLDOC01-appb-M000077
 の超球面上の一様分布を、ｐ（λ｜ａ，ｂ）は（ａ，ｂ）でパラメトライズされる逆ガンマ分布、すなわち
Figure JPOXMLDOC01-appb-M000078
 とする。ａはスパース性を調整するハイパーパラメータで、ａ＝１＋１０－１６、ｂは経験的にＡｉｍの期待値Ｅ（Ａｉｍ）から見積もれて、
Figure JPOXMLDOC01-appb-M000079
 と関係づけられる。これは更にＥ（Ｘｉｄ）と
Figure JPOXMLDOC01-appb-M000080
 なる関係があり、Ｅ（Ｘｉｄ）をＸの平均μＸで近似すると、
Figure JPOXMLDOC01-appb-M000081
 よりｂは、
Figure JPOXMLDOC01-appb-M000082
 と定まる。したがって負対数尤度関数は、
Figure JPOXMLDOC01-appb-M000083
 と書ける。これはλに対して下に凸な関数なので、
Figure JPOXMLDOC01-appb-M000084
 すなわち、
Figure JPOXMLDOC01-appb-M000085
 としてλに関する更新式が得られる。ここまでは志賀らの報告したARD-SO-NMFの導出に完全に従った。
　ノイズに対するｉ．ｉ．ｄガウス分布の仮定はＭＳデータについては当てはまらず、大きいシグナルほど大きなノイズを有する傾向があることが知られている。そこで、ノイズ分布
Figure JPOXMLDOC01-appb-M000086
 を自然同位体ピークによる線形回帰の残渣成分として求めた。すなわち、
Figure JPOXMLDOC01-appb-M000087
 ここで
Figure JPOXMLDOC01-appb-M000088
 はチャネルｄを中心にした［ｄ－３０，ｄ－２０，ｄ－１０，ｄ＋１０，ｄ＋２０，ｄ＋３０］なるチャネルで、チャネルｄの±３ｍ／ｚの同位体ピークである（チャネル間隔が０．１ｍ／ｚであることに留意）。チャネルごとの分散共分散行列
Figure JPOXMLDOC01-appb-M000089
 は
Figure JPOXMLDOC01-appb-M000090
 として求まる。これを用いて尤度関数を書き直すと、
Figure JPOXMLDOC01-appb-M000091
 となる。Rは定数行列であることに注意する。全体の負対数尤度関数は、
Figure JPOXMLDOC01-appb-M000092
  
となる。λについての更新式Ｅｑ．Ｓ３を代入し、定数項を落として簡略化すると、
Figure JPOXMLDOC01-appb-M000093
 となる。この関数のＡ、Ｓに関する最小化はｈｉｅｒａｒｃｈｉｃａｌ　ａｌｔｅｒｎａｔｉｎｇ　ｌｅａｓｔ　ｓｑｕａｒｅ　（ＨＡＬＳ）により行う。簡便のため、以下のベクトル表記を用いる。
Figure JPOXMLDOC01-appb-M000094
 ＨＡＬＳでは残差Ｘ－ＡＳを
Figure JPOXMLDOC01-appb-M000095
 として表す。ここで
Figure JPOXMLDOC01-appb-M000096
  
である。これによりＬ（Ｘ，Ａ，Ｓ，λ）は成分ｍごとに分離して書けて、
Figure JPOXMLDOC01-appb-M000097
 Ｓについてのソフト直交制約はペナルティー項
Figure JPOXMLDOC01-appb-M000098
 として目的関数Ｌに組み込むことができる。ここで
Figure JPOXMLDOC01-appb-M000099
 である。この項はｍ－成分とその他成分の非直交性を表し、ξｍは厳密に直交性が満たされる際のラグランジェ未定乗数を表す。この項は更にｗＯ∈［０，１］によって緩和される。したがって最小化すべき目的関数は、
Figure JPOXMLDOC01-appb-M000100
 となり、ａｍとｓｍに関する勾配は、
Figure JPOXMLDOC01-appb-M000101
 となり、これらをゼロと置けば、更新式は
Figure JPOXMLDOC01-appb-M000102
 となる。ｓｍは更新毎に規格化されるので、定数係数は省略した。非負制約は、更新毎に非負象限に射影することにより満たした。具体的には例えば、
Figure JPOXMLDOC01-appb-M000103
 によって射影できる。ここで
Figure JPOXMLDOC01-appb-M000104
 は要素ごとに絶対値を取ったベクトルを意味する。未定乗数ξｍは
Figure JPOXMLDOC01-appb-M000105
 をＥｑ．Ｓ６の左側から乗じ、厳密な直交条件、
Figure JPOXMLDOC01-appb-M000106
 とｗＯ＝１を適用することで、
Figure JPOXMLDOC01-appb-M000107
 と得られる。以上Ｅｑ．Ｓ１－Ｓ７を用いて以下のアルゴリズム１を提案する。
Figure JPOXMLDOC01-appb-M000108
　ここで、ＡとＳの更新毎に、似た成分同士をマージングするようにする。同様のスペクトルを有する成分同士をマージすることが志賀らに提案されているが、ここではさらに、同様の強度分布を有する成分同士もマージングすることにした。これにより、同位体ピーク、異なるイオン付加によりイオン化されたフラグメントシリーズ、モノマー数違いのオリゴマーピークシリーズなどを一つの成分とすることで、より解釈性の高い結果を与えることができる。
（Ｃａｎｏｎｉｃａｌ　ｃｏｒｒｅｌａｔｉｏｎ　ａｎａｌｙｓｉｓ　ｆｉｌｔｅｒ　（ＣＣＡ－ｆｉｌｔｅｒ）の導出）
入力：１回目ＮＭＦの出力
Figure JPOXMLDOC01-appb-M000109
 、バックグラウンドスペクトル
Figure JPOXMLDOC01-appb-M000110
出力：バックグランド由来と判定された成分番号のリスト
　１回目のＮＭＦで得られたＭ－成分の中には、バックグランドや混入物に由来する成分が含まれ、これらは組成分析結果をゆがめる可能性があるため、ＡやＳから取り除くことが好ましい。Ｍ′－成分がバックグラウンド由来の成分だとＣＣＡ－ｆｉｌｔｅｒにより判定されたとすると、
Figure JPOXMLDOC01-appb-M000111
 及び
Figure JPOXMLDOC01-appb-M000112
 はそれぞれ、
Figure JPOXMLDOC01-appb-M000113
 と
Figure JPOXMLDOC01-appb-M000114
 になる。単純のため、ＣＣＡ－ｆｉｌｔｅｒ適用後の成分数はＭ－Ｍ′になるが、一貫してＭを用いる。
　ＣＣＡ－ｆｉｌｔｅｒを用いるためには、バックグラウンドスペクトル
Figure JPOXMLDOC01-appb-M000115
 をデータセットに含めて、サンプルスペクトルと併せて１回目のＮＭＦを行う必要がある。もし何らかの混入物が予想できる場合は、その混入物について測定したスペクトルをＸＢＧとして用いることができる。ＣＣＡ－ｆｉｌｔｅｒは成分ｍ＝１，…，Ｍについてひとつずつ、そのスペクトルＳｍ：がＸＢＧに含まれているか確認していく。
　最初のステップは、Ｍ－スペクトルからなる
Figure JPOXMLDOC01-appb-M000116
 されたＳを、Ｓｍ：に似たスペクトルセット
Figure JPOXMLDOC01-appb-M000117
 と、似ていないスペクトルセット
Figure JPOXMLDOC01-appb-M000118
 に分割することである。この分割は、
Figure JPOXMLDOC01-appb-M000119
 が満たされるように行う。ｔ１∈［０，１］はある閾値で本発明では一貫してｔ１＝０．２を用いた。なお、Ｓｍ：は常にＹに含まれるので、Ｙの一列目にＳｍ：を格納した。ＺはＸＢＧと結合し、
Figure JPOXMLDOC01-appb-M000120
 とした。ＹとＺは、平均ゼロとなるように、
Figure JPOXMLDOC01-appb-M000121
 とした。これら二つのスペクトルセットに対して、ＣＣＡを行う。ＣＣＡは、二つのスペクトルセット内の線形結合により、なるべく似たペアのスペクトルを作り出す。
Figure JPOXMLDOC01-appb-M000122
及び
Figure JPOXMLDOC01-appb-M000123
 の線形結合の係数がベクトル
Figure JPOXMLDOC01-appb-M000124
 及び
Figure JPOXMLDOC01-appb-M000125
 に格納されているとする。従ってスペクトルペアはそれぞれ、
Figure JPOXMLDOC01-appb-M000126
 及び
Figure JPOXMLDOC01-appb-M000127
 と書ける。その類似度は相関係数ρによって評価され、
Figure JPOXMLDOC01-appb-M000128
 ここで、
Figure JPOXMLDOC01-appb-M000129
 　ＣＣＡの問題設定は、
Figure JPOXMLDOC01-appb-M000130
 と書ける。その解
Figure JPOXMLDOC01-appb-M000131
 は一般固有値問題の解として与えられて、
Figure JPOXMLDOC01-appb-M000132
 ここで
Figure JPOXMLDOC01-appb-M000133
 は固有ベクトル
Figure JPOXMLDOC01-appb-M000134
 をその固有値が大きい順に列ベクトルとして並べて得られる行列である。
　各固有値
Figure JPOXMLDOC01-appb-M000135
 は、対応する（ｕ＊，ｖ＊）を係数として線形結合で作られたｙとｚの相関係数になっている。ここではρ＞ｔ２を満たすすべての固有値に対応するｕ＊を取り出し、もしＳｍ：の係数に対応する１番目の要素がｕ＊に大きく寄与する大きな成分であれば、すなわち
Figure JPOXMLDOC01-appb-M000136
 であれば、成分ｍはバックグラウンドに由来する成分と判定して、システムから除く。ここでｔ２∈［０．９，０．９９］、ｔ３∈［０，１］である。これを以下にアルゴリズム２としてまとめた。
Figure JPOXMLDOC01-appb-M000137
　バックグラウンド成分が特定されたあとは、対応するＡの列ベクトル、及びＳの行ベクトルを削除することで、システムから除いた。
　ＣＣＡ－ｆｉｌｔｅｒから出力された、バックグラウンドに由来するＭ′－成分除去済みの
Figure JPOXMLDOC01-appb-M000138
 （以下、簡単のためＭ－Ｍ′をＭと置き換えて、
Figure JPOXMLDOC01-appb-M000139
 と表す）は、サンプル量及び内部標準ピークで割り付け強度を補正する。その後、サンプルｎに関するＡの小行列部分Ａ（ｎ）を一次元ベクトル化した
Figure JPOXMLDOC01-appb-M000140
 をサンプルｎの特徴ベクトルとして二回目のＮＭＦの入力としてもよいが、組成分析においてはＭ－フラグメントの温度分布は不要な情報のため、サンプル毎にすべての温度帯域を足し合わせてサンプル毎のＭ－フラグメント存在強度（ｆｒａｇｍｅｎｔ　ａｂｕｎｄａｎｃｅ；　ＦＡ）を表す
Figure JPOXMLDOC01-appb-M000141
 を入力としてもよい。ここでは簡単のため
Figure JPOXMLDOC01-appb-M000142
 について二回目のＮＭＦを行うことにする。
　以下のセクションでは、Ｎｏｎ－ｎｅｇａｔｉｖｅ　ｌｅａｓｔ　ｓｑｕａｒｅ　（ＮＮＬＳ）によるフィッティングを頻繁に使う。この問題は、最適な非負係数
Figure JPOXMLDOC01-appb-M000143
 を用いて、定数行列
Figure JPOXMLDOC01-appb-M000144
 の列ベクトル線形結合により
Figure JPOXMLDOC01-appb-M000145
 を近似するもので、
Figure JPOXMLDOC01-appb-M000146
 により求められる。これを解くには、ａｌｔｅｒｎａｔｉｎｇ　ｄｉｒｅｃｔｉｏｎ　ｍｕｌｔｉｐｌｉｅｒ　ｍｅｔｈｏｄｓ　（ＡＤＭＭ）（１）を含む数多くの最適化手法が使えるが、ここではＦｕ．　ｅｔ．　ａｌ．　（２）により開発されたＡＤＭＭ－ＮＮＬＳを用いてこの問題をとき、解を
Figure JPOXMLDOC01-appb-M000147
 と表記する。近似するベクトルセットＹ＝［ｙ１，．．，ｙＬ］に対応する非負係数ベクトルｘｌ ＊（ｌ＝１，．．，Ｌ）は個別に計算できて、
Figure JPOXMLDOC01-appb-M000148
 または行列形式で、
Figure JPOXMLDOC01-appb-M000149
 と書くことにする。ここで
Figure JPOXMLDOC01-appb-M000150
 である。同様の問題で、係数ベクトルの総和が１になるような拘束条件付の問題はｆｕｌｌｙ　ｃｏｎｓｔｒａｉｎｅｄ　ｌｅａｓｔ　ｓｑｕａｒｅ　（ＦＣＬＳ）とよばれ、ＡＤＭＭを用いて同じく解くことができて、
Figure JPOXMLDOC01-appb-M000151
 または行列表記で、
Figure JPOXMLDOC01-appb-M000152
 と書くことにする。
（２回目のＮＭＦ）
入力：サンプル毎のＦＡを表す
Figure JPOXMLDOC01-appb-M000153
 、基底Ｍ－フラグメントごとのスペクトル
Figure JPOXMLDOC01-appb-M000154
 、ポリマー成分数K
出力：サンプル毎のポリマー重量分率
Figure JPOXMLDOC01-appb-M000155
 、単位重量Ｋ－基底ポリマーごとのＦＡ
Figure JPOXMLDOC01-appb-M000156
　２回目のＮＭＦ：
Figure JPOXMLDOC01-appb-M000157
 は、リーマン計量を用いてその近似残差
Figure JPOXMLDOC01-appb-M000158
 を評価する、すなわち、
Figure JPOXMLDOC01-appb-M000159
 ここで、
Figure JPOXMLDOC01-appb-M000160
 である。下三角行列
Figure JPOXMLDOC01-appb-M000161
 をコレスキー分解Ｇ＝ＬＬＴにより得れば、
Figure JPOXMLDOC01-appb-M000162
 と書ける。ここで
Figure JPOXMLDOC01-appb-M000163
 であり、
Figure JPOXMLDOC01-appb-M000164
 と行列分解することと等値である。
Figure JPOXMLDOC01-appb-M000165
 の行ベクトルで張られるシンプレックスの体積項を
Figure JPOXMLDOC01-appb-M000166
 の行ベクトル間の非直交性項を
Figure JPOXMLDOC01-appb-M000167
 として、α＞０とβ∈［０，１］を重みとすれば、
Figure JPOXMLDOC01-appb-M000168
 により
Figure JPOXMLDOC01-appb-M000169
 が求まる。ここで、Ｆｕらは（３）外れ値へのロバスト性を付与するために、
Figure JPOXMLDOC01-appb-M000170
 なる重み行列の導入を提案している。ここで、ｐ∈（０，２］、かつεは小さな正則化パラメータである。本発明では、一貫してｐ＝１．５かつε＝１０－８とした。最適化問題は
Figure JPOXMLDOC01-appb-M000171
 となる。Ｃと
Figure JPOXMLDOC01-appb-M000172
 はｂｌｏｃｋ　ｃｏｏｒｄｉｎａｔｅ　ｄｅｓｃｅｎｔ　（ＢＣＤ）　理論により交互に最適化する。今、ｔ回の最適化を行い、
Figure JPOXMLDOC01-appb-M000173
 を得たとする。なお、初期化には、ｖｅｒｔｅｘ　ｃｏｍｐｏｎｅｎｔ　ａｎａｌｙｓｉｓ　（ＶＣＡ）（４）を用いて、
Figure JPOXMLDOC01-appb-M000174
 のＮ－行ベクトルの中からシンプレックスのＫ－頂点に近いＫ個の行ベクトルを選択し、これを
Figure JPOXMLDOC01-appb-M000175
 とした。
Figure JPOXMLDOC01-appb-M000176
 に基づくＣのアップデートは、
Figure JPOXMLDOC01-appb-M000177
 により単純に行える。そこでＣ（ｔ）に基づく
Figure JPOXMLDOC01-appb-M000178
 から
Figure JPOXMLDOC01-appb-M000179
 へのアップデートについて考える。なお、
Figure JPOXMLDOC01-appb-M000180
 で、
Figure JPOXMLDOC01-appb-M000181
 は未定乗数である。また
Figure JPOXMLDOC01-appb-M000182
 ここで、
Figure JPOXMLDOC01-appb-M000183
 は小さな正則化パラメータで、一貫してτ＝１０－８を用いた。ここで、
Figure JPOXMLDOC01-appb-M000184
 のＭａｊｏｒｉｚｅｒ関数を導入する。接線不等式より、一期前の
Figure JPOXMLDOC01-appb-M000185
 に基づいて、
Figure JPOXMLDOC01-appb-M000186
 が成立する。ここで
Figure JPOXMLDOC01-appb-M000187
 は
Figure JPOXMLDOC01-appb-M000188
 のＨ（ｔ）における勾配で、
Figure JPOXMLDOC01-appb-M000189
 より、
Figure JPOXMLDOC01-appb-M000190
 ここで、
Figure JPOXMLDOC01-appb-M000191
 かつ、ｃｏｎｓｔは定数項である。したがって、
Figure JPOXMLDOC01-appb-M000192
 を
Figure JPOXMLDOC01-appb-M000193
 で置き換えて、すべてのペナルティー項を合わせて書くと、
Figure JPOXMLDOC01-appb-M000194
 ここで
Figure JPOXMLDOC01-appb-M000195
 従って
Figure JPOXMLDOC01-appb-M000196
 のアップデートは以下を解くことで得られる。
Figure JPOXMLDOC01-appb-M000197
 複雑な制約条件を整理するために、
Figure JPOXMLDOC01-appb-M000198
 なる制約をラグランジェ未定乗数
Figure JPOXMLDOC01-appb-M000199
 を用いて、目的関数に組み込み、ＡＤＭＭの枠組みでこの問題を解く。
Figure JPOXMLDOC01-appb-M000200
 ここで、
Figure JPOXMLDOC01-appb-M000201
 かつμはＡＤＭＭのハイパーパラメータである（ここでは一貫してμ＝１とした）。Ｚ′は一期前のＺを表し、
Figure JPOXMLDOC01-appb-M000202
 のときに目的関数はＺについて最大化された、
Figure JPOXMLDOC01-appb-M000203
 を与える。これを
Figure JPOXMLDOC01-appb-M000204
 について最適化すればよい。
Figure JPOXMLDOC01-appb-M000205
 をサイクリックに最適化し、アルゴリズムをAlgorithm3にまとめた。
Figure JPOXMLDOC01-appb-M000206
目的関数
Figure JPOXMLDOC01-appb-M000207
 は
Figure JPOXMLDOC01-appb-M000208
 については拘束条件なしの二次関数の最適化なので、
Figure JPOXMLDOC01-appb-M000209
 より最小化される。また
Figure JPOXMLDOC01-appb-M000210
 はＮＮＬＳにより解けて、
Figure JPOXMLDOC01-appb-M000211
 と更新される。以上より元の問題Ｅｑ．Ｓ１１を解くための
Figure JPOXMLDOC01-appb-M000212
 の更新式が得られた。最後にＶ＝αＦ（ｔ）＋βΛ（ｔ）の更新について考える。Ｆ（ｔ）については
Figure JPOXMLDOC01-appb-M000213
 と簡単に求まる。Λ（ｔ）は厳密な直交条件、すなわち
Figure JPOXMLDOC01-appb-M000214
 と
Figure JPOXMLDOC01-appb-M000215
 を組み合わせて、
Figure JPOXMLDOC01-appb-M000216
 をＥｑ．Ｓ１９の右側からかけると、
Figure JPOXMLDOC01-appb-M000217
 と得られる。Ｅｑ．Ｓ１８とＥｑ．Ｓ２０より、
Figure JPOXMLDOC01-appb-M000218
 が与えられる。以上を用いて、問題Ｅｑ．Ｓ１１を解くためのアルゴリズムは以下のとおりである。
Figure JPOXMLDOC01-appb-M000219
（ＳとＢの張る超平面へのテストデータの射影）
　データセットからＳとＢを推定したのちに、ＳとＢの推定には用いなかったデータ（ここではテストデータと呼ぶ）
Figure JPOXMLDOC01-appb-M000220
 を、ＳとＢの張る超平面へ射影する方法について述べる。
　まずＳ－超平面への射影により、
Figure JPOXMLDOC01-appb-M000221
 が得られる。すなわち、
Figure JPOXMLDOC01-appb-M000222
 温度帯域について総和をとり、
Figure JPOXMLDOC01-appb-M000223
 に変換したのちに、Ｂ－超平面への射影と規格化により、
Figure JPOXMLDOC01-appb-M000224
 が得られる。すなわち、
Figure JPOXMLDOC01-appb-M000225
　なお、データセットとハイパーパラメータは以下のとおりである。
Figure JPOXMLDOC01-appb-T000226
　また、引用文献は、以下のとおりである。
Figure JPOXMLDOC01-appb-T000227
（組成推定装置）
　次に、本発明の実施形態に係る組成推定装置について、図面を参照しながら説明する。図３は、本発明の実施形態に係る組成推定装置のハードウェア構成図である。
　組成推定装置３０は、質量分析装置３１と、情報処理装置３３とを備え、質量分析装置３１は、減圧ポンプ３２を備えている。情報処理装置３３は、プロセッサ３４と、記憶デバイス３５と、表示デバイス３６と、入力デバイス３７とを有し、質量分析装置３１と情報処理装置３３とは相互にデータを送受信できるよう構成されている。
　プロセッサ３４は、例えば、マイクロプロセッサ、プロセッサコア、マルチプロセッサ、ＡＳＩＣ（ａｐｐｌｉｃａｔｉｏｎ－ｓｐｅｃｉｆｉｃ　ｉｎｔｅｇｒａｔｅｄ　ｃｉｒｃｕｉｔ）、ＦＰＧＡ（ｆｉｅｌｄ　ｐｒｏｇｒａｍｍａｂｌｅ　ｇａｔｅ　ａｒｒａｙ）、及び、ＧＰＧＰＵ（Ｇｅｎｅｒａｌ－ｐｕｒｐｏｓｅ　ｃｏｍｐｕｔｉｎｇ　ｏｎ　ｇｒａｐｈｉｃｓ　ｐｒｏｃｅｓｓｉｎｇ　ｕｎｉｔｓ）等である。
　記憶デバイス３５は、各種プログラム、及び、データを一時的に、及び／又は、非一時的に記憶する機能を有し、プロセッサ３４の作業エリアを提供する。
　記憶デバイス３５は、例えば、ＲＯＭ（Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）、ＲＡＭ（Ｒａｎｄｏｍ　Ａｃｃｅｓｓ　Ｍｅｍｏｒｙ）、ＨＤＤ（Ｈａｒｄ　Ｄｉｓｋ　Ｄｒｉｖｅ）、フラッシュメモリ、及び、ＳＳＤ（Ｓｏｌｉｄ　Ｓｔａｔｅ　Ｄｒｉｖｅ）等である。
　入力デバイス３７は、各種情報入力を受け付け、また、組成推定装置３０への指示の入力を受け付けることができる。入力デバイス３７は、キーボード、マウス、スキャナ、及び、タッチパネル等でよい。
　表示デバイス３６は、組成推定装置３０のステータス、分析の進捗、及び、組成推定結果等を表示できる。表示デバイス３６は、液晶ディスプレイ、及び、有機ＥＬ（Ｅｌｅｃｔｒｏ　Ｌｕｍｉｎｅｓｃｅｎｃｅ）ディスプレイ等でよい。
　また、表示デバイス３６は、入力デバイス３７と一体として構成されていてもよい。この場合、表示デバイス３６がタッチパネルディスプレイであって、ＧＵＩ（Ｇｒａｐｈｉｃａｌ　Ｕｓｅｒ　Ｉｎｔｅｒｆａｃｅ）を提供する形態が挙げられる。
　データバスにより相互にデータを通信可能なプロセッサ３４、記憶デバイス３５、入力デバイス３７、及び、表示デバイス３６を備える、情報処理装置３３は、典型的にはコンピュータである。
　質量分析装置３１は、減圧ポンプ３２を備え、典型的には、ＤＡＲＴイオン源、検体加熱装置、及び、飛行時間型質量分析計を備える質量分析装置である。質量分析装置３１が備えるイオン源、及び、質量分析計は、いずれも非限定的な例であり、本組成推定装置が備える質量分析装置の構成は上記に限定されない。
　図４は、本発明の実施形態に係る組成推定装置３０の機能ブロック図である。組成推定装置３０は、全検体の質量分析を行う質量分析部４１と、質量分析部４１によって得られたマススペクトルを処理する情報処理部４２とを有する。
　質量分析部４１は、質量分析装置３１を含んで構成され、質量分析装置３１が後述する制御部４３に制御されて実現する機能である。
　質量分析部４１は、ロードされた検体セット（検体、学習用検体、及び、バックグラウンド検体、図４中、符号ａ）について、加熱しながら熱脱着、及び／又は、熱分解されて生ずるガス成分をイオン化し、順次質量分析して、マススペクトルを出力する（図４中、符号ｃ）。
　この際、質量分析装置３１の設置されている環境に起因した環境要因（図４中、符号ｂ、例えば、減圧ポンプ３２から生ずる（ポリ）ジメチルシロキサン）も同時に質量分析装置３１にロードされ、マススペクトルに影響を与える。
　次に、情報処理装置３３を含んで構成される情報処理部４２の各部の機能について説明する。
　制御部４３は、プロセッサ３４を含んで構成される。制御部４３は、各部を制御して、組成推定装置３０の各機能を実現する。
　記憶部４４は、記憶デバイス３５を含んで構成される。記憶部４４には予め各部の制御のためのプログラムが記憶され、また、各部によって入出力されるデータの一時的、又は、非一時的な記憶領域を提供する。記憶部４４は、制御部４３によって制御され、データを読み出したり、書き込んだりすることができる。
　入力部４５は入力デバイス３７を含んで構成される。また、出力部４６は、表示デバイス３６を含んで構成される。制御部４３がこれらを制御することで、組成推定装置３０の使用者からの入力を受け付けて、これを記憶部４４に記憶させたり、記憶部４４に記憶された情報を各部に送信したり、出力部４６を介して組成推定装置３０の使用者に対して表示したりできる。
　データ行列作成部４７は、記憶部４４に記憶されたプログラムが制御部４３により実行されて実現される機能である。
　データ行列作成部４７は、質量分析部４１によって加熱温度ごとに取得されたマススペクトル（図４中、符号ｃ）を各行に格納した二次元マススペクトルからデータ行列Ｘを作成し、後述するＮＭＦ処理部４８に渡す（図４中、符号ｄ）。なお、データ行列Ｘの詳細については、すでに説明したとおりである。
　ＮＭＦ処理部４８は、記憶部４４に記憶されたプログラムが制御部４３により実行されて実現される機能である。
　ＮＭＦ処理部４８は、データ行列作成部４７によって作成されたデータ行列Ｘ（図４中、符号ｄ）を非負値行列因子分解（ＮＭＦ）し、規格化された基底スペクトル行列（Ｓ）とその強度分布行列（Ｃ）との積（Ｘ＝Ｃ×ＳＴ）に分解する。その結果は、補正処理部４９に渡される（図４中、符号ｅ）。なお、ＮＭＦ処理の詳細はすでに説明したとおりである。
　補正処理部４９は、記憶部４４に記憶されたプログラムが制御部４３により実行されて実現される機能である。
　補正処理部４９は、ＮＭＦ処理部４８によって作成された基底スペクトル行列とデータ行列との正準相関分析によって強度分布行列を補正して、補正後強度分布行列を作成する機能である。作成された補正後強度分布行列はベクトル処理部５０に渡される（図４中、符号ｆ）。
　なお、補正処理部４９は、上記の正準相関分析による補正に加えて、すでに説明した強度補正を更に行ってもよい。この際、強度補正に用いるキャリブレーション検体の組成（ｃａｌｉｂ検体組成）がある場合は、入力部４５を介して補正処理部４９に渡される（図４中、符号ｋ）。正準相関分析による補正（ＣＣＡ　ｆｉｌｔｅｒ）、及び、強度補正の詳細はすでに説明したとおりである。
　ベクトル処理部５０は、記憶部４４に記憶されたプログラムが制御部４３により実行されて実現される機能である。
　ベクトル処理部５０は、補正処理部４９により作成された補正後強度分布行列を全検体のそれぞれに対応する小行列に分割し、これを特徴ベクトルとしてベクトル空間内で全検体を表現する。この結果は、エンドメンバー決定部５１に渡される（図４中、符号ｇ）。
　エンドメンバー決定部５１は、記憶部４４に記憶されたプログラムが制御部４３により実行されて実現される機能である。
　エンドメンバー決定部５１は、ベクトル処理部５０によって作成された特徴ベクトルの全てを内包するＫ－１次元単体を設定し、Ｋ－１次元単体におけるＫ個のエンドメンバーを決定する。この結果は、含有量比計算部５２に渡される（図４中、符号ｈ）。
　この際、学習用検体にエンドメンバーが含まれる場合であって、判別ラベルを有する場合には、入力部４５から判別ラベルの情報がエンドメンバー決定部５１に渡される（図４中、符号ｊ）。
　しかし、すでに説明したとおり、エンドメンバーの決定（エンドメンバー決定部５１の処理）には判別ラベルは必須ではない。
　含有量比計算部５２は、記憶部４４に記憶されたプログラムが制御部４３により実行されて実現される機能である。
　含有量比計算部５２は、エンドメンバー決定部５１によって決定されたＫ個のエンドメンバーと、ベクトル処理部５０により作成された推定対象検体の特徴ベクトルとのユークリッド距離をそれぞれ計算し、Ｋ個のエンドメンバーと検体の特徴ベクトルとのユークリッド距離の比から、推定対象検体中のＫ個の成分の含有量比を推定する。この含有量比の推定結果（図４中、符号ｉ）は、出力部４６から出力される。
　なお、この際、学習用検体の組成推定結果もあわせて出力してもよい。
　本組成推定装置よれば、溶媒に溶解しにくく、従来の組成分析が行いにくい固体の検体であっても、特段の前処理なしに分析に供することができ、かつ、正確な定量分析結果が得られる。本組成推定装置は、特に、複数の成分の混合物である検体の品質管理、及び、不良原因の調査等に適用した場合、結果を得るために必要な時間を大幅に短縮することができる。
　以下では、実施例をもとに、本組成分析方法について詳述する。なお、本発明の組成分析方法は下記の特定の実施例によって限定的に解釈されるものでない。
［実施例１］
（検体の準備）
　ポリメチルメタクリレート（ｐｏｌｙ（ｍｅｔｈｙｌ　ｍｅｔｈａｃｒｙｌａｔｅ）、記号「Ｍ」）、ポリスチレン（ｐｏｌｙｓｔｙｒｅｎｅ、記号「Ｓ」）、ポリエチルメタクリレート（ｐｏｌｙ（ｅｔｈｙｌ　ｍｅｔｈａｃｒｙｌａｔｅ）、記号「Ｅ」）の三種類の高分子のジオキサン溶液（約４質量％、高分子の固形分として０．４ｍｇ程度）を、サンプルポット上にキャストした後に一晩風乾して検体を得た。表４は、各検体の組成である。なお、質量分析は、ＳＨＩＭＡＤＺＵ社製ＬＣＭＳ－２０２０を検出器として、ＩｏｎＳｅｎｓｅ社製ＤＡＲＴイオン源として、バイオクロマト社製イオンロケットをサンプルの加熱部として用いて行った。
Figure JPOXMLDOC01-appb-T000228
　得られたマススペクトルについて、以下の（ｉ）～（ｉｖ）を実装したアルゴリズムによって解析を行った。
（ｉ）Ｃ行列ＡＲＤとＳ行列へのソフト直交制約
（ｉｉ）ＮＭＦ最適化中のマージング条件
（ｉｉｉ）ＣＣＡ－ｆｉｌｔｅｒによるＮＭＦアーティファクト、及び、バックグラウンド、コンタミネーション除去
（ｉｖ）環境変数を補正項にいれることで、キャリブレーション不要の定量分析
　解析に使用したハイパーパラメータは以下のとおりである。
ＮＭＦ：
ｏｒｔｈｏｇｏｎａｌ　ｃｏｎｓｔｒａｉｎ　ｗｅｉｇｈｔ　ｗ＝０．２，　ｉｎｉｔｉａｌ　ｃｏｍｐｏｎｅｎｔ　ｎｕｍｂｅｒ　Ｋ＝５０，　ｓｍａｌｌ　ｐｏｓｉｔｉｖｅ　ｎｕｍｂｅｒ　ａ＝１＋０．１－１５，　ｉｔｅｒａｔｉｏｎ　ｎｕｍｂｅｒ　Ｌ＝２０００，　ｔｈｒｅｓｈｏｌｄ＝０．９
ＣＣＡ－Ｆｉｌｔｅｒ：
ｔ１＝０．８，ｔ２＝０．２，
　図５は、「ＭＳＥ＿ｃａｌｉｂ」についてのオリジナルなスペクトル（Ａ）と、上記処理を行った後、再構築したスペクトル（Ｂ）を表す図である。両者の波形は類似していることがわかる。すなわち、上記処理によってオリジナルなスペクトルの特徴は保存され、組成分析に用いやすいデータとなったことが示された。
　また、オリジナルデータにはエアコン等から排出されるフタル酸エステル等のバックグラウンドが３９０ｍ／ｚに観測される（図５（Ａ）中、「ＢＧ」と表示した部分）が、再構築されたピーク（Ｂ）では、これが除去されている。また、検体の作成に用いられた溶媒であるジオキサン（図５（Ａ）中、「ＤＯＸ」と表示した部分）も除去されていた。これらはいずれもノイズに該当するものであり、本方法によって、効果的にノイズが除去されたことが示された。
　また、図６は、ＮＭＦ分解後のＣ行列の一部Ｃδ（ＭＳＥ＿Ｃａｌｉｂ）と、ＣＣＡ－Ｆｉｌｔｅｒによってバックグラウンド及び不純物のスペクトルを除いた後のＳＴである。
　図７は、Ｍ、Ｓ、Ｅ三成分系の組成解析結果である。図７中、中空円で表され、凡例に「Ｃａｌｉｂｒａｔｅｄ」とあるのは、ＭＳＥ＿ｃａｌｉｂをキャリブレーション検体とした場合、つまり、ＭＳＥ＿ｃａｌｉｂの組成情報を既知として与え、強度分布行列を強度補正した場合の検体の組成推定結果である。
　また、図７中、中実円で表され、凡例に「Ｎｏｎ－ｃａｌｉｂｒａｔｅｄ」とあるのは、環境変数（すなわち、６１０ｍ／ｚ：ポリジメチルシロキサンのピーク強度）と検体質量の積を用いて、強度分布行列を強度補正した場合の検体の組成推定結果である。
　また、図７中、星印で表され、凡例に「Ｋｎｏｗｎ」とあるのは、実際の組成（真の値）を表している。
　図７に示した結果から、「Ｃａｌｉｂｒａｔｅｄ」、「Ｎｏｎ－ｃａｌｉｂｒａｔｅｄ」のいずれの場合にも真の値に近い、優れた推定結果が得られた。
　次に、図８は、ＮＭＦの条件から、（１）の直交制約を解除した場合の組成推定結果である。各プロットの定義は、図７と同様である。
　図８の結果から、ＮＭＦが直交制約を有する方が、より優れた推定結果が得られることがわかった。
　次に、図９は、ＭＮＦ最適化の条件から、ｍｅｒｇｉｎｇ　Ｓのみ解除した場合の組成推定結果である。各プロットの定義は、図７と同様である。
　図９の結果から、ＭＮＦ最適化の条件がｍｅｒｇｉｎｇ　Ｓを有する方が、より優れた推定結果が得られることがわかった。
　なお、成分数はＡＲＤによって３６までしか減らせないが、Ｍｅｒｇｉｎｇ　Ｓの条件を維持したときは成分数９でより精度よく組成が推定できることから、Ｍｅｒｇｉｎｇ　Ｓの重要性がわかった。
　図１０は、ＣＣＡ－ｆｉｌｔｅｒを解除した場合の結果（比較例）の組成推定結果である。各プロットの定義は図７と同様である。
　図１０においては、検体「ＳＥ」は、Ｍの含有量は０％にも関わらず１％程度入っているという計算結果となった。０％を１％と推定してしまうと、数％の推定精度もないということになり実用的ではない。
　図１１は、「ＭＳＥ＿ｃａｌｉｂ」についてのオリジナルなスペクトル（上、ｒａｗ　ｏｆ　ＭＳＥ＿ｃａｌｉｂ）と、ＣＣＡ－ｆｉｌｔｅｒを解除した上記で処理を行った後、再構築したスペクトル（下、ｒｅｃｏｎｓｔｒｕｃｔｅｄ　ｏｆ　ＭＳＥ＿ｃａｌｉｂ）を表す図である。
　この結果から、ＣＣＡ－ｆｉｌｔｅｒを解除すると、ＤＯＸ、ＢＧ（図５の説明を参照）のいずれのノイズも除去できていないことがわかった。
　図１２は、環境変数による補正を行わなかった場合の組成推定結果である。図１２の結果においては、中実円で表された「Ｎｏｎ－ｃａｌｉｂｒａｔｅｄ」がＫｎｏｗｎからのずれが大きく成っており、環境変数による補正により、得られる推定結果がより正確となることがわかった。
［実施例２］
（エンドメンバーを含まない検体による推定）
　学習用検体を含む検体（群）が、エンドメンバーを含まない場合にも精度よく推定できることを確かめるために、表５に記載したサンプルを調製して実施例１と同様の試験を行った。
　なお、表中、記号「Ｍ」「Ｓ」「Ｅ」が意味するものは、表４（実施例１）と同様である。また、各サンプルは、単独の成分（ポリマー）の重量分率が８０％を超えないように調整されている。
Figure JPOXMLDOC01-appb-T000229
　上記検体について、実施例１と同様の方法で二次元マススペクトルを得て、データ行列に変換した。得られたデータ行列について、すでに説明した１回目のＮＭＦ処理、すなわち、
Figure JPOXMLDOC01-appb-M000230
 を実施した。なお、ノイズの分散共分散行列の見積もりは、Ｘのみを用いて計算できるので、１回目のＮＭＦ処理におけるモジュール内部処理とした。
　次に、Ａに対するＣＣＡ－Ｆｉｌｔｅｒ、重さや標準ピークによる強度補正を行った。なお、組成分析には温度分布情報は不要なため、温度軸に沿ってサンプル毎に総和を取って、
Figure JPOXMLDOC01-appb-M000231
 としてもよい。
　次に、
Figure JPOXMLDOC01-appb-M000232
 の２回目のＮＭＦ処理を実施した。なお、リーマン計量による評価等は、２回目のＮＭＦ処理の内部処理として実施した。図１５はその結果として、ＣとＢＳとを図示したものである。
　図中Ｃでは、本方法によって、推定された組成を丸印で、使用したサンプルの組成（正しい組成）を星印で示している。いずれもよい一致を示している。
　ここでは、使用した検体にエンドメンバー、すなわち、Ｅ、Ｓ、Ｍのピュアポリマーが含まれていないことが、Ｃの各頂点にプロットがないことからもわかる。
　また、図中ＢＳは、推測によって得られたＥ、Ｍ、Ｓのピュアポリマーのスペクトルと、実測のＥ、Ｍ、Ｓのスペクトルとを重ねて示しているが、両者に差はなく、絶対強度を含めて、よく推定されていることが確かめられた。
　このように、本発明の方法によれば、図１５のＢＳに示されたように、Ｅ、Ｍ、Ｓのピュアポリマーのスペクトル、及び、混合物の組成が、正確に推測され得ることが確かめられた。
　なお、テストデータがある場合には、ＸｔｅｓｔをＳとＢに連続射影してＣｔｅｓｔを得る。
［実施例３］
（エンドメンバーを１つ含む検体による推定）
　次に、（学習用）検体のうち１つがエンドメンバーを含む場合に、他の検体における成分含有量が任意となる（Ｋ－１次元単体の内接円の外側の領域でなくてもよい）ことを検証した。
　表５に記載したサンプルのうち、Ｓを４０質量％以上含むサンプルを除き、かつ、サンプルとしてＭピュアポリマーを加えたこと以外は、実施例２と同様にして、組成推定を行った。図１６は結果である。
　図１６において、中実円で表された点は、Ｂ、Ｓの推定に用いた検体の推定結果を表し、星印がその実際の組成である。この結果から、Ｓについては、内接円の外側の領域に検体はなかったものの（中空円は、推定に使用しなかった検体である）、精度よく推定することができることが確認された。
　上記原因のひとつは、すでに説明した２回目のＮＭＦ処理（Ｅｑ．１１）において基底の非直交性項
Figure JPOXMLDOC01-appb-M000233
 をペナルティー項として組み込んだことによるものと考えらえる。
　材料開発の多くの場合で、材料性能を各構成成分に還元して考察することが求められる。これには、成分の同定と定量が必要だが、系が複雑化するほどその難易度は上がる。異なる成分のピークが互いにオーバーラップする可能性が低い高解像度の質量分析は、複雑な混合物系における成分の化学構造を同定するうえで最も強力な分析手法と言える。一方で、質量分析は定量性を有さないことが問題となる。また、本質的には分離を必要としないものの、その結果得られるスペクトルは複雑で現実的には読解不能に陥りやすい。
　本発明では、材料の破片を大気下・前処理なしで直接熱分解させ、そこから生じるガスをＤＡＲＴ－ＭＳで分析し、そのスペクトルを教師無し学習することで、各成分のピュアなエンドメンバースペクトルを推定し、その線形和で検体のスペクトルを表す手法を開発した。これは、複雑な組成からなる材料であっても分離を必要とせずに、その成分の同定と定量を可能とする。
　同定については、他の成分に由来するピークをすべて排除した可読性の高いエンドメンバースペクトルを得られれば、あとは比較的容易である。一方で、定量分析の妨げとなる、ピークごとにイオン化効率が異なりピーク強度とその存在比が比例しない、という事実は本発明においても同様である。しかし、本発明ではエンドメンバースペクトルの各ピーク強度を、規格化された相対値ではなく、絶対強度として推定できるため、エンドメンバースペクトルの線形和で検体スペクトルを表現する際、その線形係数をそのまま組成比として求めることができる。
　検体を成分ごとに分解し、その含有比まで求める本発明は、材料の性能を発揮する要因となったキー成分の同定など、材料開発において重宝されるツールになり得るのは明らかであるが、さらには劣化成分の含有量検査等の品質管理、危険物含有量等の法令順守検査、及び、権利侵害の検出のための検査等、きわめて広範な応用が期待される。さらには樹脂の劣化・寿命予測、樹脂・有機材料における極微不純物成分の検出、嗅覚や味覚の定量分析、食品等のトレーサビリティー・産地偽装防止タグなどの技術が加速すると期待できる。
　また、Ｐｌａｓｔｉｃｓ　Ｅｕｒｏｐｅは２０２１年９月に欧州委員会に対して２０３０年までにプラスチック容器包装に再生プラスチック３０％含有の義務化を提言するなど、再生材の利用が促進されているが、含有量を測定する一般的手法は確立されていない。本発明はバージン剤と再生材の比率を正確に測定し得る技術としても、品質管理や法令順守検査などへの応用が期待される。
１０：Ｋ－１次元単体、　１２：内接超球、　１３－１５：エンドメンバー、　１６－１８：学習用検体、　１９：領域、　３０：組成推定装置、　３１：質量分析装置　、３２：減圧ポンプ、　３３：情報処理装置、　３４：プロセッサ、　３５：記憶デバイス、　３６：表示デバイス、　３７：入力デバイス、　４１：質量分析部、　４２：情報処理部、　４３：制御部、　４４：記憶部、　４５：入力部、　４６：出力部、　４７：データ行列作成部、　４８：ＮＭＦ処理部、　４９：補正処理部、　５０：ベクトル処理部、　５１：エンドメンバー決定部、　５２：含有量比計算部

 

Claims (21)

  1. 　Ｋを１以上の整数としたとき、Ｋ種類の成分から選択される少なくとも１種の前記成分を含む推定対象検体における、前記成分の含有量比を推定する方法であって、
    　前記Ｋ種類の成分から選択される少なくとも１種の成分を含み、互いに組成の異なるＫ個以上の学習用検体、及び、前記成分を含まないバックグラウンド検体を準備することと、
    　前記推定対象検体、前記学習用検体、及び、バックグラウンド検体を含む検体セットに含まれるそれぞれの検体を加熱しながら、熱脱着、及び／又は、熱分解によって生ずるガス成分を順次イオン化し、マススペクトルを連続的に観測することと、
    　加熱温度ごとに得られた前記マススペクトルを各行に格納し、前記検体ごとの二次元マススペクトルを得て、前記二次元マススペクトルの少なくとも２つ以上をまとめてデータ行列に変換することと、
    　前記データ行列を非負値行列因子分解し、規格化された基底スペクトル行列とその強度分布行列の積に分解するＮＭＦ処理を行うことと、
    　前記基底スペクトル行列と前記データ行列との正準相関分析によって強度分布行列に含まれるノイズ成分を抽出し、前記ノイズ成分の影響を減ずるよう前記強度分布行列を補正して、補正後強度分布行列を得ることと、
    　前記補正後強度分布行列を前記検体のそれぞれに対応する小行列に分割し、前記小行列を特徴ベクトルとして、ベクトル空間内で前記検体のそれぞれを表現することと、
    　前記特徴ベクトルの全てを内包するＫ－１次元単体を設定し、前記Ｋ－１次元単体におけるＫ個のエンドメンバーを決定することと、
    　前記Ｋ個のエンドメンバーと、前記推定対象検体の特徴ベクトルとの、ユークリッド距離をそれぞれ計算し、前記ユークリッド距離の比により、前記推定対象検体中の前記成分の含有量比を推定することと、を含み、
    　前記Ｋが３以上の場合、前記Ｋ－１次元単体に内接する超球体の外側の領域のそれぞれに、前記学習用検体の特徴ベクトルの少なくとも１つが位置する、又は、前記学習用検体が前記エンドメンバーの少なくとも１つを含む、方法。
  2. 　前記Ｋが２以上の整数であり、前記推定対象検体が、前記成分の混合物である、請求項１に記載の方法。
  3. 　前記学習用検体が、前記エンドメンバーを含む、請求項１又は２に記載の方法。
  4. 　前記エンドメンバーの決定が、前記学習用検体に付された判別ラベルに基づき行われる、請求項３に記載の方法。
  5. 　前記エンドメンバーの決定が、前記学習用検体の前記特徴ベクトルの頂点成分分析によって行われる、請求項３に記載の方法。
  6. 　前記エンドメンバーの決定が、前記Ｋ－１次元単体が最小体積となるように頂点を定めるアルゴリズムによって実施される、請求項１に記載の方法。
  7. 　前記エンドメンバーの決定が、前記補正後強度分布行列を非負値行列因子分解し、前記検体中のＫ種類の前記成分の重量分率を表す行列と、Ｋ種類の個別の前記成分ごとのフラグメント存在強度を表す行列との積に分解する第２のＮＭＦ処理により行われる、請求項６に記載の方法。
  8. 　前記学習用検体が、前記エンドメンバーを含まない、請求項６又は７に記載の方法。
  9. 　前記補正後強度分布行列を得ることは、更に、前記強度分布行列の強度補正を行うことを含む、請求項１～８のいずれか１項に記載の方法。
  10. 　前記検体セットが、更に、キャリブレーション検体を含み、
    　前記キャリブレーション検体は、前記Ｋ種類の成分を全て含み、かつ、その組成が既知であり、
    　前記強度補正は、前記強度分布行列を規格化することを含む、請求項９に記載の方法。
  11. 　前記強度補正は、前記強度分布行列の少なくとも一部を、前記検体セットのうち対応する前記検体の質量と環境変数の積で割り付けることを含み、
    　前記環境変数は、前記観測の際における前記成分のイオン化効率への影響を表す変数である、請求項９に記載の方法。
  12. 　前記環境変数が、前記観測の際に、雰囲気中に所定量含まれる分子量５０～１５００の化合物、又は、前記検体セットのそれぞれに所定量含まれる分子量５０～５００の有機低分子化合物のマススペクトルのピークの合計値である、請求項１１に記載の方法。
  13. 　Ｋを１以上の整数としたとき、Ｋ種類の成分から選択される少なくとも１種の前記成分を含む推定対象検体における、前記成分の含有量比を推定する組成推定装置であって、
    　前記Ｋ種類の成分から選択される少なくとも１種の成分を含み、互いに組成の異なるＫ個以上の学習用検体、前記成分を含まないバックグラウンド検体、及び、前記推定対象検体を含む検体セットに含まれるそれぞれの検体を加熱しながら、熱脱着、及び／又は、熱分解によって生ずるガス成分を順次イオン化し、マススペクトルを連続的に観測する質量分析装置と、
    　前記観測されたマススペクトルを処理する情報処理装置と、を備え、
    　前記情報処理装置は、
    　加熱温度ごとに得られた前記マススペクトルを各行に格納し、前記検体ごとの二次元マススペクトルを得て、前記二次元マススペクトルの少なくとも２つ以上をまとめてデータ行列に変換する、データ行列作成部と、
    　前記データ行列を非負値行列因子分解し、規格化された基底スペクトル行列とその強度分布行列の積に分解するＮＭＦ処理を行う、ＮＭＦ処理部と、
    　前記基底スペクトル行列と前記データ行列との正準相関分析によって強度分布行列に含まれるノイズ成分を抽出し、前記ノイズ成分の影響を減ずるよう前記強度分布行列を補正して、補正後強度分布行列を作成する、補正処理部と、
    　前記補正後強度分布行列を前記検体セットのそれぞれに対応する小行列に分割し、前記小行列を特徴ベクトルとして、ベクトル空間内で前記検体のそれぞれを表現する、ベクトル処理部と、
    　前記特徴ベクトルの全てを内包するＫ－１次元単体を設定し、前記Ｋ－１次元単体におけるＫ個のエンドメンバーを決定する、エンドメンバー決定部と、
    　前記Ｋ個のエンドメンバーと、前記推定対象検体の特徴ベクトルとの、ユークリッド距離をそれぞれ計算し、前記ユークリッド距離の比により、前記推定対象検体中の前記成分の含有量比を推定する、含有量比計算部と、を含み、
    　前記Ｋが３以上の場合、前記Ｋ－１次元単体に内接する超球体の外側の領域のそれぞれに、前記学習用検体の特徴ベクトルの少なくとも１つが位置する、又は、前記学習用検体が前記エンドメンバーの少なくとも１つを含む、組成推定装置。
  14. 　前記エンドメンバー決定部は、前記学習用検体に付された判別ラベルに基づいて前記エンドメンバーを決定する、請求項１３に記載の組成推定装置。
  15. 　前記エンドメンバー決定部は、前記Ｋ－１次元単体が最小体積となるように頂点を定めるアルゴリズムによって前記エンドメンバーを決定する、請求項１３に記載の組成推定装置。
  16. 　前記エンドメンバー決定部は、前記補正後強度分布行列を非負値行列因子分解し、前記検体中のＫ種類の前記成分の重量分率を表す行列と、Ｋ種類の個別の前記成分ごとのフラグメント存在強度を表す行列との積に分解する第２のＮＭＦ処理により前記エンドメンバーを決定する、請求項１５に記載の組成推定装置。
  17. 　前記補正処理部は更に、前記強度分布行列の強度補正を行う、請求項１５に記載の組成推定装置。
  18. 　前記検体セットが、更に、キャリブレーション検体を含み、
    　前記キャリブレーション検体が、前記Ｋ種類の成分を全て含み、かつ、その組成が既知であり、
    　前記強度補正が、前記強度分布行列を規格化することを含む、請求項１７に記載の組成推定装置。
  19. 　前記強度補正は、前記強度分布行列の少なくとも一部を、前記検体セットのうち対応する前記検体の質量と環境変数の積で割り付けることを含み、
    　前記環境変数は、前記観測の際における前記成分のイオン化効率への影響を表す変数である、請求項１７に記載の組成推定装置。
  20. 　前記質量分析装置は減圧ポンプを更に備え、
    　前記環境変数が、前記減圧ポンプの稼働によって生ずる物質のマススペクトルのピークの合計値である、請求項１９に記載の組成推定装置。
  21. 　Ｋを１以上の整数としたとき、Ｋ種類の成分から選択される少なくとも１種の前記成分を含む推定対象検体における、前記成分の含有量比を推定する組成推定装置において用いられるプログラムであって、
    　前記Ｋ種類の成分から選択される少なくとも１種の成分を含み、互いに組成の異なるＫ個以上の学習用検体、前記成分を含まないバックグラウンド検体、及び、前記推定対象検体を含む検体セットに含まれるそれぞれの検体を加熱しながら、熱脱着、及び／又は、熱分解によって生ずるガス成分を順次イオン化し、マススペクトルを連続的に観測する質量分析装置と、
    　前記観測されたマススペクトルを処理する情報処理装置と、を備え、
    　前記質量分析装置によって、加熱温度ごとに取得された前記マススペクトルを各行に格納し、前記検体ごとの二次元マススペクトルを得て、前記二次元マススペクトルの少なくとも２つ以上をまとめてデータ行列に変換する、データ行列作成機能と、
    　前記データ行列を非負値行列因子分解し、規格化された基底スペクトル行列とその強度分布行列の積に分解するＮＭＦ処理を行う、ＮＭＦ処理機能と、
    　前記基底スペクトル行列と前記データ行列との正準相関分析によって強度分布行列に含まれるノイズ成分を抽出し、前記ノイズ成分の影響を減ずるよう前記強度分布行列を補正して、補正後強度分布行列を作成する、補正処理機能と、
    　前記補正後強度分布行列を前記検体のそれぞれに対応する小行列に分割し、前記小行列を特徴ベクトルとして、ベクトル空間内で前記検体のそれぞれを表現する、ベクトル処理機能と、
    　前記特徴ベクトルの全てを内包するＫ－１次元単体を設定し、前記Ｋ－１次元単体におけるＫ個のエンドメンバーを決定する、エンドメンバー決定機能と、
    　前記Ｋ個のエンドメンバーと、前記推定対象検体の特徴ベクトルとの、ユークリッド距離をそれぞれ計算し、前記ユークリッド距離の比により、前記推定対象検体中の前記成分の含有量比を推定する、含有量比計算機能と、を含み、
    　前記Ｋが３以上の場合、前記Ｋ－１次元単体に内接する超球体の外側の領域のそれぞれに、前記学習用検体の特徴ベクトルの少なくとも１つが位置する、又は、前記学習用検体が前記エンドメンバーの少なくとも１つを含む、プログラム。

     
PCT/JP2022/022813 2021-06-24 2022-06-06 検体に含まれる成分の含有量比の推定方法、組成推定装置、及び、プログラム Ceased WO2022270289A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2023529793A JP7628731B2 (ja) 2021-06-24 2022-06-06 検体に含まれる成分の含有量比の推定方法、組成推定装置、及び、プログラム
US18/569,560 US20240297031A1 (en) 2021-06-24 2022-06-06 Method of inferring content ratio of component in sample, composition inference device, and program
EP22828206.7A EP4361624B1 (en) 2021-06-24 2022-06-06 Method for estimating content ratio of components contained in sample, composition estimating device, and program

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2021104510 2021-06-24
JP2021-104510 2021-06-24

Publications (1)

Publication Number Publication Date
WO2022270289A1 true WO2022270289A1 (ja) 2022-12-29

Family

ID=84544540

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2022/022813 Ceased WO2022270289A1 (ja) 2021-06-24 2022-06-06 検体に含まれる成分の含有量比の推定方法、組成推定装置、及び、プログラム

Country Status (4)

Country Link
US (1) US20240297031A1 (ja)
EP (1) EP4361624B1 (ja)
JP (1) JP7628731B2 (ja)
WO (1) WO2022270289A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2025033010A1 (ja) 2023-08-09 2025-02-13 国立研究開発法人物質・材料研究機構 組成推定方法、組成推定装置、プログラム、検体容器、及び、熱重量・質量分析方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012516445A (ja) * 2009-01-27 2012-07-19 ホロジック，インコーポレイテッド 体液における新生児敗血症の検出のためのバイオマーカー
JP2013528287A (ja) * 2010-06-10 2013-07-08 インターナショナル・ビジネス・マシーンズ・コーポレーション 質量スペクトルを分析するための方法、コンピュータ・プログラム、およびシステム
WO2016120958A1 (ja) * 2015-01-26 2016-08-04 株式会社島津製作所 ３次元スペクトルデータ処理装置及び処理方法
WO2018158801A1 (ja) * 2017-02-28 2018-09-07 株式会社島津製作所 スペクトルデータの特徴抽出装置および方法
JP2021081365A (ja) * 2019-11-21 2021-05-27 株式会社島津製作所 糖ペプチド解析装置

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1997049B1 (en) * 2006-03-21 2019-03-13 Metabolon Inc. A system, method, and computer program product for analyzing spectrometry data to indentify and quantify individual components in a sample
WO2007112597A1 (en) * 2006-04-05 2007-10-11 National Research Council Of Canada Blind extraction of pure component mass spectra from overlapping mass spectrometric peaks

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012516445A (ja) * 2009-01-27 2012-07-19 ホロジック，インコーポレイテッド 体液における新生児敗血症の検出のためのバイオマーカー
JP2013528287A (ja) * 2010-06-10 2013-07-08 インターナショナル・ビジネス・マシーンズ・コーポレーション 質量スペクトルを分析するための方法、コンピュータ・プログラム、およびシステム
WO2016120958A1 (ja) * 2015-01-26 2016-08-04 株式会社島津製作所 ３次元スペクトルデータ処理装置及び処理方法
WO2018158801A1 (ja) * 2017-02-28 2018-09-07 株式会社島津製作所 スペクトルデータの特徴抽出装置および方法
JP2021081365A (ja) * 2019-11-21 2021-05-27 株式会社島津製作所 糖ペプチド解析装置

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
"Robust Volume Minimization-based Matrix Factorization for remote Sensing and Document Clustering", IEEE TRANSACTIONS ON SIGNAL PROCESSING, vol. 64, no. 23, 1 December 2016 (2016-12-01)
ANALYTICAL CHEMISTRY, vol. 92, 2020, pages 1925 - 1933
J. M. B. DIASJ. M. P. NASCIMENTO: "Vertex component analysis: A geometric-based approach to unmix hyperspectral data", SIGNAL IMAGE PROCESS. REMOTE SENS., vol. 43, 2006, pages 415 - 439
S. BOYDN. PARIKHE. CHUB. PELEATOJ. ECKSTEIN: "Distributed optimization and statistical learning via the alternating direction method of multipliers", FOUND. TRENDS MACH. LEARN., vol. 3, 2010, pages 1 - 122, XP055127725, DOI: 10.1561/2200000016
See also references of EP4361624A4
SHIGA, M. ET AL.: "Sparse modeling of EELS and EDX spectral imaging data by nonnegative matrix factorization", ULTRAMICROSCOPY, vol. 170, 2016, pages 43 - 59
X. FUK. HUANGB. YANGW. K. MAN. D. SIDIROPOULOS: "Robust Volume Minimization-Based Matrix Factorization for Remote Sensing and Document Clustering", IEEE TRANS. SIGNAL PROCESS., vol. 64, 2016, pages 6254 - 6268, XP011625333, DOI: 10.1109/TSP.2016.2602800

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2025033010A1 (ja) 2023-08-09 2025-02-13 国立研究開発法人物質・材料研究機構 組成推定方法、組成推定装置、プログラム、検体容器、及び、熱重量・質量分析方法

Also Published As

Publication number Publication date
EP4361624A1 (en) 2024-05-01
JPWO2022270289A1 (ja) 2022-12-29
EP4361624A4 (en) 2025-07-02
EP4361624B1 (en) 2026-05-06
US20240297031A1 (en) 2024-09-05
JP7628731B2 (ja) 2025-02-12

Similar Documents

Publication Publication Date Title
Burger et al. Data handling in hyperspectral image analysis
JP6089345B2 (ja) 時および／または空間系列ファイルの多成分回帰／多成分分析
Gorrochategui et al. Chemometric strategy for untargeted lipidomics: biomarker detection and identification in stressed human placental cells
Trindade et al. Non-negative matrix factorisation of large mass spectrometry datasets
CN103959426B (zh) 用于通过质谱术识别微生物的方法
Moorthy et al. Pattern similarity measures applied to mass spectra
CN109738413B (zh) 基于稀疏非负最小二乘的混合物拉曼光谱定性分析方法
JP6748085B2 (ja) 干渉検出および着目ピークのデコンボルーション
Lewis et al. Comparing multiple statistical methods for inverse prediction in nuclear forensics applications
Koo et al. Analysis of metabolomic profiling data acquired on GC–MS
Borsoi et al. Fast unmixing and change detection in multitemporal hyperspectral data
Franceschi et al. Self‐organizing maps: A versatile tool for the automatic analysis of untargeted imaging datasets
US20240393310A1 (en) Biomarker reflectance signatures for disease detection and classification
Mani-Varnosfaderani et al. Estimating complicated baselines in analytical signals using the iterative training of Bayesian regularized artificial neural networks
Gibb et al. Mass spectrometry analysis using MALDIquant
Rapin et al. Application of non-negative matrix factorization to LC/MS data
Lalor et al. Sim-to-real supervised domain adaptation for radioisotope identification
Deng et al. An infinite Gaussian mixture model with its application in hyperspectral unmixing
JP7757246B2 (ja) 試料分析装置及び方法
Teurtrie et al. From STEM-EDXS data to phase separation and quantification using physics-guided NMF
EP4361624B1 (en) Method for estimating content ratio of components contained in sample, composition estimating device, and program
Zhang et al. Chromatographic fingerprinting and chemometric techniques for quality control of herb medicines
JP2008500537A (ja) 分光計によって生成されたデータからスペクトルを抽出するためのシステム及び方法
Negoita et al. Artificial intelligence application designed to screen for new psychoactive drugs based on their ATR-FTIR spectra
Caiafa On the conditions for valid objective functions in blind separation of independent and dependent sources

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 22828206

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 2023529793

Country of ref document: JP

WWE Wipo information: entry into national phase

Ref document number: 18569560

Country of ref document: US

WWE Wipo information: entry into national phase

Ref document number: 2022828206

Country of ref document: EP

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 2022828206

Country of ref document: EP

Effective date: 20240124

WWG Wipo information: grant in national office

Ref document number: 2022828206

Country of ref document: EP