org.apache.commons.math.stat.regression.SimpleRegression Maven / Gradle / Ivy
Show all versions of commons-math Show documentation
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.stat.regression;
import java.io.Serializable;
import org.apache.commons.math.MathException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.distribution.TDistribution;
import org.apache.commons.math.distribution.TDistributionImpl;
/**
* Estimates an ordinary least squares regression model
* with one independent variable.
*
* y = intercept + slope * x
*
* Standard errors for intercept and slope are
* available as well as ANOVA, r-square and Pearson's r statistics.
*
* Observations (x,y pairs) can be added to the model one at a time or they
* can be provided in a 2-dimensional array. The observations are not stored
* in memory, so there is no limit to the number of observations that can be
* added to the model.
*
* Usage Notes:
* - When there are fewer than two observations in the model, or when
* there is no variation in the x values (i.e. all x values are the same)
* all statistics return
NaN. At least two observations with
* different x coordinates are requred to estimate a bivariate regression
* model.
*
* - getters for the statistics always compute values based on the current
* set of observations -- i.e., you can get statistics, then add more data
* and get updated statistics without using a new instance. There is no
* "compute" method that updates all statistics. Each of the getters performs
* the necessary computations to return the requested statistic.
*
*
* @version $Revision: 773189 $ $Date: 2009-05-09 05:57:04 -0400 (Sat, 09 May 2009) $
*/
public class SimpleRegression implements Serializable {
/** Serializable version identifier */
private static final long serialVersionUID = -3004689053607543335L;
/** the distribution used to compute inference statistics. */
private TDistribution distribution;
/** sum of x values */
private double sumX = 0d;
/** total variation in x (sum of squared deviations from xbar) */
private double sumXX = 0d;
/** sum of y values */
private double sumY = 0d;
/** total variation in y (sum of squared deviations from ybar) */
private double sumYY = 0d;
/** sum of products */
private double sumXY = 0d;
/** number of observations */
private long n = 0;
/** mean of accumulated x values, used in updating formulas */
private double xbar = 0;
/** mean of accumulated y values, used in updating formulas */
private double ybar = 0;
// ---------------------Public methods--------------------------------------
/**
* Create an empty SimpleRegression instance
*/
public SimpleRegression() {
this(new TDistributionImpl(1.0));
}
/**
* Create an empty SimpleRegression using the given distribution object to
* compute inference statistics.
* @param t the distribution used to compute inference statistics.
* @since 1.2
*/
public SimpleRegression(TDistribution t) {
super();
setDistribution(t);
}
/**
* Adds the observation (x,y) to the regression data set.
*
* Uses updating formulas for means and sums of squares defined in
* "Algorithms for Computing the Sample Variance: Analysis and
* Recommendations", Chan, T.F., Golub, G.H., and LeVeque, R.J.
* 1983, American Statistician, vol. 37, pp. 242-247, referenced in
* Weisberg, S. "Applied Linear Regression". 2nd Ed. 1985.
*
*
* @param x independent variable value
* @param y dependent variable value
*/
public void addData(double x, double y) {
if (n == 0) {
xbar = x;
ybar = y;
} else {
double dx = x - xbar;
double dy = y - ybar;
sumXX += dx * dx * (double) n / (n + 1d);
sumYY += dy * dy * (double) n / (n + 1d);
sumXY += dx * dy * (double) n / (n + 1d);
xbar += dx / (n + 1.0);
ybar += dy / (n + 1.0);
}
sumX += x;
sumY += y;
n++;
if (n > 2) {
distribution.setDegreesOfFreedom(n - 2);
}
}
/**
* Removes the observation (x,y) from the regression data set.
*
* Mirrors the addData method. This method permits the use of
* SimpleRegression instances in streaming mode where the regression
* is applied to a sliding "window" of observations, however the caller is
* responsible for maintaining the set of observations in the window.
*
* The method has no effect if there are no points of data (i.e. n=0)
*
* @param x independent variable value
* @param y dependent variable value
*/
public void removeData(double x, double y) {
if (n > 0) {
double dx = x - xbar;
double dy = y - ybar;
sumXX -= dx * dx * (double) n / (n - 1d);
sumYY -= dy * dy * (double) n / (n - 1d);
sumXY -= dx * dy * (double) n / (n - 1d);
xbar -= dx / (n - 1.0);
ybar -= dy / (n - 1.0);
sumX -= x;
sumY -= y;
n--;
if (n > 2) {
distribution.setDegreesOfFreedom(n - 2);
}
}
}
/**
* Adds the observations represented by the elements in
* data.
*
* (data[0][0],data[0][1]) will be the first observation, then
* (data[1][0],data[1][1]), etc.
*
* This method does not replace data that has already been added. The
* observations represented by data are added to the existing
* dataset.
*
* To replace all data, use clear() before adding the new
* data.
*
* @param data array of observations to be added
*/
public void addData(double[][] data) {
for (int i = 0; i < data.length; i++) {
addData(data[i][0], data[i][1]);
}
}
/**
* Removes observations represented by the elements in data.
*
* If the array is larger than the current n, only the first n elements are
* processed. This method permits the use of SimpleRegression instances in
* streaming mode where the regression is applied to a sliding "window" of
* observations, however the caller is responsible for maintaining the set
* of observations in the window.
*
* To remove all data, use clear().
*
* @param data array of observations to be removed
*/
public void removeData(double[][] data) {
for (int i = 0; i < data.length && n > 0; i++) {
removeData(data[i][0], data[i][1]);
}
}
/**
* Clears all data from the model.
*/
public void clear() {
sumX = 0d;
sumXX = 0d;
sumY = 0d;
sumYY = 0d;
sumXY = 0d;
n = 0;
}
/**
* Returns the number of observations that have been added to the model.
*
* @return n number of observations that have been added.
*/
public long getN() {
return n;
}
/**
* Returns the "predicted" y value associated with the
* supplied x value, based on the data that has been
* added to the model when this method is activated.
*
* predict(x) = intercept + slope * x
*
* Preconditions:
* - At least two observations (with at least two different x values)
* must have been added before invoking this method. If this method is
* invoked before a model can be estimated,
Double,NaN is
* returned.
*
*
* @param x input x value
* @return predicted y value
*/
public double predict(double x) {
double b1 = getSlope();
return getIntercept(b1) + b1 * x;
}
/**
* Returns the intercept of the estimated regression line.
*
* The least squares estimate of the intercept is computed using the
* normal equations.
* The intercept is sometimes denoted b0.
*
* Preconditions:
* - At least two observations (with at least two different x values)
* must have been added before invoking this method. If this method is
* invoked before a model can be estimated,
Double,NaN is
* returned.
*
*
* @return the intercept of the regression line
*/
public double getIntercept() {
return getIntercept(getSlope());
}
/**
* Returns the slope of the estimated regression line.
*
* The least squares estimate of the slope is computed using the
* normal equations.
* The slope is sometimes denoted b1.
*
* Preconditions:
* - At least two observations (with at least two different x values)
* must have been added before invoking this method. If this method is
* invoked before a model can be estimated,
Double.NaN is
* returned.
*
*
* @return the slope of the regression line
*/
public double getSlope() {
if (n < 2) {
return Double.NaN; //not enough data
}
if (Math.abs(sumXX) < 10 * Double.MIN_VALUE) {
return Double.NaN; //not enough variation in x
}
return sumXY / sumXX;
}
/**
* Returns the
* sum of squared errors (SSE) associated with the regression
* model.
*
* The sum is computed using the computational formula
*
* SSE = SYY - (SXY * SXY / SXX)
*
* where SYY is the sum of the squared deviations of the y
* values about their mean, SXX is similarly defined and
* SXY is the sum of the products of x and y mean deviations.
*
* The sums are accumulated using the updating algorithm referenced in
* {@link #addData}.
*
* The return value is constrained to be non-negative - i.e., if due to
* rounding errors the computational formula returns a negative result,
* 0 is returned.
*
* Preconditions:
* - At least two observations (with at least two different x values)
* must have been added before invoking this method. If this method is
* invoked before a model can be estimated,
Double,NaN is
* returned.
*
*
* @return sum of squared errors associated with the regression model
*/
public double getSumSquaredErrors() {
return Math.max(0d, sumYY - sumXY * sumXY / sumXX);
}
/**
* Returns the sum of squared deviations of the y values about their mean.
*
* This is defined as SSTO
* here.
*
* If n < 2, this returns Double.NaN.
*
* @return sum of squared deviations of y values
*/
public double getTotalSumSquares() {
if (n < 2) {
return Double.NaN;
}
return sumYY;
}
/**
* Returns the sum of squared deviations of the x values about their mean.
*
* If n < 2, this returns Double.NaN.
*
* @return sum of squared deviations of x values
*/
public double getXSumSquares() {
if (n < 2) {
return Double.NaN;
}
return sumXX;
}
/**
* Returns the sum of crossproducts, xi*yi.
*
* @return sum of cross products
*/
public double getSumOfCrossProducts() {
return sumXY;
}
/**
* Returns the sum of squared deviations of the predicted y values about
* their mean (which equals the mean of y).
*
* This is usually abbreviated SSR or SSM. It is defined as SSM
* here
*
* Preconditions:
* - At least two observations (with at least two different x values)
* must have been added before invoking this method. If this method is
* invoked before a model can be estimated,
Double.NaN is
* returned.
*
*
* @return sum of squared deviations of predicted y values
*/
public double getRegressionSumSquares() {
return getRegressionSumSquares(getSlope());
}
/**
* Returns the sum of squared errors divided by the degrees of freedom,
* usually abbreviated MSE.
*
* If there are fewer than three data pairs in the model,
* or if there is no variation in x, this returns
* Double.NaN.
*
* @return sum of squared deviations of y values
*/
public double getMeanSquareError() {
if (n < 3) {
return Double.NaN;
}
return getSumSquaredErrors() / (n - 2);
}
/**
* Returns
* Pearson's product moment correlation coefficient,
* usually denoted r.
*
* Preconditions:
* - At least two observations (with at least two different x values)
* must have been added before invoking this method. If this method is
* invoked before a model can be estimated,
Double,NaN is
* returned.
*
*
* @return Pearson's r
*/
public double getR() {
double b1 = getSlope();
double result = Math.sqrt(getRSquare());
if (b1 < 0) {
result = -result;
}
return result;
}
/**
* Returns the
* coefficient of determination,
* usually denoted r-square.
*
* Preconditions:
* - At least two observations (with at least two different x values)
* must have been added before invoking this method. If this method is
* invoked before a model can be estimated,
Double,NaN is
* returned.
*
*
* @return r-square
*/
public double getRSquare() {
double ssto = getTotalSumSquares();
return (ssto - getSumSquaredErrors()) / ssto;
}
/**
* Returns the
* standard error of the intercept estimate,
* usually denoted s(b0).
*
* If there are fewer that three observations in the
* model, or if there is no variation in x, this returns
* Double.NaN.
*
* @return standard error associated with intercept estimate
*/
public double getInterceptStdErr() {
return Math.sqrt(
getMeanSquareError() * ((1d / (double) n) + (xbar * xbar) / sumXX));
}
/**
* Returns the standard
* error of the slope estimate,
* usually denoted s(b1).
*
* If there are fewer that three data pairs in the model,
* or if there is no variation in x, this returns Double.NaN.
*
*
* @return standard error associated with slope estimate
*/
public double getSlopeStdErr() {
return Math.sqrt(getMeanSquareError() / sumXX);
}
/**
* Returns the half-width of a 95% confidence interval for the slope
* estimate.
*
* The 95% confidence interval is
*
* (getSlope() - getSlopeConfidenceInterval(),
* getSlope() + getSlopeConfidenceInterval())
*
* If there are fewer that three observations in the
* model, or if there is no variation in x, this returns
* Double.NaN.
*
* Usage Note:
* The validity of this statistic depends on the assumption that the
* observations included in the model are drawn from a
*
* Bivariate Normal Distribution.
*
* @return half-width of 95% confidence interval for the slope estimate
* @throws MathException if the confidence interval can not be computed.
*/
public double getSlopeConfidenceInterval() throws MathException {
return getSlopeConfidenceInterval(0.05d);
}
/**
* Returns the half-width of a (100-100*alpha)% confidence interval for
* the slope estimate.
*
* The (100-100*alpha)% confidence interval is
*
* (getSlope() - getSlopeConfidenceInterval(),
* getSlope() + getSlopeConfidenceInterval())
*
* To request, for example, a 99% confidence interval, use
* alpha = .01
*
* Usage Note:
* The validity of this statistic depends on the assumption that the
* observations included in the model are drawn from a
*
* Bivariate Normal Distribution.
*
* Preconditions:
* - If there are fewer that three observations in the
* model, or if there is no variation in x, this returns
*
Double.NaN.
*
* (0 < alpha < 1); otherwise an
* IllegalArgumentException is thrown.
*
*
* @param alpha the desired significance level
* @return half-width of 95% confidence interval for the slope estimate
* @throws MathException if the confidence interval can not be computed.
*/
public double getSlopeConfidenceInterval(double alpha)
throws MathException {
if (alpha >= 1 || alpha <= 0) {
throw MathRuntimeException.createIllegalArgumentException(
"out of bounds significance level {0}, must be between {1} and {2}",
alpha, 0.0, 1.0);
}
return getSlopeStdErr() *
distribution.inverseCumulativeProbability(1d - alpha / 2d);
}
/**
* Returns the significance level of the slope (equiv) correlation.
*
* Specifically, the returned value is the smallest alpha
* such that the slope confidence interval with significance level
* equal to alpha does not include 0.
* On regression output, this is often denoted Prob(|t| > 0)
*
* Usage Note:
* The validity of this statistic depends on the assumption that the
* observations included in the model are drawn from a
*
* Bivariate Normal Distribution.
*
* If there are fewer that three observations in the
* model, or if there is no variation in x, this returns
* Double.NaN.
*
* @return significance level for slope/correlation
* @throws MathException if the significance level can not be computed.
*/
public double getSignificance() throws MathException {
return 2d * (1.0 - distribution.cumulativeProbability(
Math.abs(getSlope()) / getSlopeStdErr()));
}
// ---------------------Private methods-----------------------------------
/**
* Returns the intercept of the estimated regression line, given the slope.
*
* Will return NaN if slope is NaN.
*
* @param slope current slope
* @return the intercept of the regression line
*/
private double getIntercept(double slope) {
return (sumY - slope * sumX) / n;
}
/**
* Computes SSR from b1.
*
* @param slope regression slope estimate
* @return sum of squared deviations of predicted y values
*/
private double getRegressionSumSquares(double slope) {
return slope * slope * sumXX;
}
/**
* Modify the distribution used to compute inference statistics.
* @param value the new distribution
* @since 1.2
*/
public void setDistribution(TDistribution value) {
distribution = value;
// modify degrees of freedom
if (n > 2) {
distribution.setDegreesOfFreedom(n - 2);
}
}
}