lsqsn

From DavinciWiki
Jump to: navigation, search

Description

Constrained least squares fitting, where some components are allowed to be negative and others are not

Arguments and Return Values

Arguments: A two-dimensional matrix, a column vector, a scalar specifying the number of forced components, and an optional fourth input specifying the format to return

Return Value: A vector with the fitted coefficients, or a structure containing the fitted coefficients and other fit quantities

Usage

Syntax: lsqsn(library_matrix, vector_to_be_fitted, n_unc[, mode])

library_matrix - A two-dimensional matrix (dimensions [n, m, 1]) with the library of known samples (e.g. spectra from mineral library), and a column vector with the unknown sample (e.g. spectrum of some rock) to be fitted with the library.

vector_to_be_fitted - A column vector (dimensions [1, m, 1]) with the sample to be fitted.

The first two inputs must have the same y-dimension.


n_unc - A nonnegative integer specifying the number of components which are unconstrained (allowed but not required to be negative). The function assumes that the nonnegative samples are the first columns in library_matrix, and the unconstrained samples are at the end. So, for a given input of n_unc, samples 1 though (n - n_unc) must be nonnegative, and samples (n - n_unc + 1) through n are unconstrained.

This function basically sets up and runs calls to lsqnn, using binbyte to iterate through all combinations of positive and negative for the unconstrained components. Therefore increasing n_unc by one doubles the required number of computations. In order to limit computational time, n_unc must be less than or equal to 30, and values close to this limit will probably take too long to be feasible.


If using this function for spectral deconvolution, the first input should be a matrix where each column corresponds to a different known sample (usually a known mineral) and each row gives the emissivity or reflectivity value of the samples at a certain wavelength or wavenumber. The second input should be a vector with the unknown sample to be analyzed, and the returned coefficients will be the concentration of each known sample.

The fitted coefficients calculated by this function must be positive or zero for the first (n - nunc) samples; the coefficients for samples (n - nunc + 1) through n are unconstrained. This is for applications such as spectral deconvolution where it does not make physical sense to have negative coefficients for most samples, but a few, such as the blackbody term or the atmospheric components, might be negative. For cases where any of the coefficients are allowed to be negative, use lsq. For cases where none of the coefficients are allowed to be negative, use lsqnn.


An optional fourth input (passed by value) controls the output mode. If $4=0, the function will only return the fitted coefficients. If $4 is nonzero, or if it is omitted, the function will display a summary of the fit errors and will return a structure with the fitted coefficients, the fitted spectrum, the fit errors, and the number of positive endmembers after each iteration (field np).

If $4=0, lsqsn returns a column vector containing the fitted coefficients. Otherwise, it returns a structure with the fields A (fitted coefficients), Yfit (fitted spectrum), and other fields quantifying the errors in the fitted spectrum.


This function returns A such that X*A best approximates y, such that the sum of the squared errors is minimized. The data are fit to y = a1*x1 + a2*x2 + ... + an*xn, subject to the constraint that all ai≥0 for 1 ≤ i ≤ n - nunc. Note that no constant offset (e.g. blackbody for emissivity spectra) is automatically included by this function; that must be added elsewhere if needed. Many spectral fitting functions (e.g. sma) have their own option to add the constant term.


If lsqsn() is entered with no arguments, it prints its description, as shown below.

Examples

dv> lsqsn()

Constrained least squares, where some elements are allowed to be
 negative and others are not
lsqsn(X,y,n) returns A such that X*A best approximates y,
 where the last n elements are allowed to be negative but
 the first elements must be positive (or zero)
n must be less than or equal to 30
In essence, this function just sets up a bunch of calls to lsqnn,
 which effectively means that each additional unconstrained
 term effectively doubles the number of computations.
Optional fourth input controls output mode
 If $4 is 1 (or omitted), will print summary and return structure
  with errors and fitted spectrum, in addition to calculated
  coefficients.
 If $4 is 0, will only return structure with coefficients
 This means that the format of the output depends on $4.
S.Marshall 11-09-2010

dv> Xb = cat(X, clone(1., 1,923,1), axis=x) # Adding blackbody
74x923x1 array of float, bsq format [273,208 bytes]
dv> y
1x923x1 array of float, bsq format [3,692 bytes]
dv> A = lsqsn(Xb, y, 1, 0)
1x74x1 array of double, bsq format [592 bytes]
dv> A = lsqsn(Xb, y, 1, 1)
Errors in fitted spectrum:
Mean absolute error: 0.011729
RMS error: 0.019572
Maximum absolute error: 0.182021
struct, 5 elements
    A: 1x74x1 array of double, bsq format [592 bytes]
    Yfit: 1x923x1 array of double, bsq format [7,384 bytes]
    MAE: 0.0117294792933117
    RMSE: 0.0195723106839733
    maxE: 0.182020887732506
dv> A = lsqsn(Xb, y, 1)
Errors in fitted spectrum:
Mean absolute error: 0.011729
RMS error: 0.019572
Maximum absolute error: 0.182021
struct, 5 elements
    A: 1x74x1 array of double, bsq format [592 bytes]
    Yfit: 1x923x1 array of double, bsq format [7,384 bytes]
    MAE: 0.0117294792933117
    RMSE: 0.0195723106839733
    maxE: 0.182020887732506

DavinciWiki Mini-Nav Bar

Contents


Contact Developers

  • davinci-dev [AT] mars.asu.edu

All other topics

  • See navigation on the left

Related Functions

Recent Library Changes

Created On: 12-02-2010
Modified On: 12-05-2010

Personal tools