Loading...
Searching...
No Matches
EST_ols.cc
/*************************************************************************/
/* */
/* Centre for Speech Technology Research */
/* University of Edinburgh, UK */
/* Copyright (c) 1998 */
/* All Rights Reserved. */
/* */
/* Permission is hereby granted, free of charge, to use and distribute */
/* this software and its documentation without restriction, including */
/* without limitation the rights to use, copy, modify, merge, publish, */
/* distribute, sublicense, and/or sell copies of this work, and to */
/* permit persons to whom this work is furnished to do so, subject to */
/* the following conditions: */
/* 1. The code must retain the above copyright notice, this list of */
/* conditions and the following disclaimer. */
/* 2. Any modifications must be clearly marked as such. */
/* 3. Original authors' names are not deleted. */
/* 4. The authors' names are not used to endorse or promote products */
/* derived from this software without specific prior written */
/* permission. */
/* */
/* THE UNIVERSITY OF EDINBURGH AND THE CONTRIBUTORS TO THIS WORK */
/* DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING */
/* ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT */
/* SHALL THE UNIVERSITY OF EDINBURGH NOR THE CONTRIBUTORS BE LIABLE */
/* FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES */
/* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN */
/* AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, */
/* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF */
/* THIS SOFTWARE. */
/* */
/*************************************************************************/
/* Author : Alan W Black (and lots of books) */
/* Date : January 1998 */
/*-----------------------------------------------------------------------*/
/* Ordinary Least Squares/Linear regression */
/* */
/*=======================================================================*/
#include <cmath>
#include "EST_multistats.h"
#include "EST_simplestats.h"
static void ols_load_selected_feats(const EST_FMatrix &X,
const EST_IVector &included,
EST_FMatrix &Xl);
static int ols_stepwise_find_best(const EST_FMatrix &X,
const EST_FMatrix &Y,
EST_IVector &included,
EST_FMatrix &coeffs,
float &bscore,
int &best_feat,
const EST_FMatrix &Xtest,
const EST_FMatrix &Ytest,
const EST_StrList &feat_names
);
int ols(const EST_FMatrix &X,const EST_FMatrix &Y, EST_FMatrix &coeffs)
{
// Ordinary least squares, X contains the samples with 1 (for intercept)
// in column 0, Y contains the single values.
EST_FMatrix Xplus;
if (!pseudo_inverse(X,Xplus))
return FALSE;
multiply(Xplus,Y,coeffs);
return TRUE;
}
int robust_ols(const EST_FMatrix &X,
const EST_FMatrix &Y,
EST_FMatrix &coeffs)
{
EST_IVector included;
int i;
included.resize(X.num_columns());
for (i=0; i<included.length(); i++)
included.a_no_check(i) = TRUE;
return robust_ols(X,Y,included,coeffs);
}
int robust_ols(const EST_FMatrix &X,
const EST_FMatrix &Y,
EST_IVector &included,
EST_FMatrix &coeffs)
{
// as with ols but if the pseudo inverse fails remove the offending
// features and try again until it works, this can be costly but
// its saves *you* from finding the singularity
// This expands the output and puts weights of 0 for omitted features
EST_FMatrix Xl;
EST_FMatrix coeffsl;
EST_FMatrix Xplus;
int i,j,singularity=-1;
if (X.num_rows() <= X.num_columns())
{
cerr << "OLS: less rows than columns, so cannot find solution."
<< endl;
return FALSE;
}
if (X.num_columns() != included.length())
{
cerr << "OLS: `included' list wrong size: internal error."
<< endl;
return FALSE;
}
while (TRUE)
{
ols_load_selected_feats(X,included,Xl);
if (pseudo_inverse(Xl,Xplus,singularity))
{
multiply(Xplus,Y,coeffsl);
break;
}
else
{ // found a singularity so try again without that column
// remap singularity position back to X
int s;
for (s=i=0; i<singularity; i++)
{
s++;
while ((included(s) == FALSE) ||
(included(s) == OLS_IGNORE))
s++;
}
if (included(s) == FALSE)
{ // oops
cerr << "OLS: found singularity twice, shouldn't happen"
<< endl;
return FALSE;
}
else
{
cerr << "OLS: omitting singularity in column " << s << endl;
included[s] = FALSE;
}
}
}
// Map coefficients back, making coefficient 0 for singular cols
coeffs.resize(X.num_columns(),1);
for (j=i=0; i<X.num_columns(); i++)
if (included(i))
{
coeffs.a_no_check(i,0) = coeffsl(j,0);
j++;
}
else
coeffs.a_no_check(i,0) = 0.0;
return TRUE;
}
static void ols_load_selected_feats(const EST_FMatrix &X,
const EST_IVector &included,
EST_FMatrix &Xl)
{
int i,j,k,width;
for (width=i=0; i<included.length(); i++)
if (included(i) == TRUE)
width++;
Xl.resize(X.num_rows(),width);
for (i=0; i<X.num_rows(); i++)
for (k=j=0; j < X.num_columns(); j++)
if (included(j) == TRUE)
{
Xl.a_no_check(i,k) = X.a_no_check(i,j);
k++;
}
}
int ols_apply(const EST_FMatrix &samples,
const EST_FMatrix &coeffs,
EST_FMatrix &res)
{
// Apply coefficients to samples for res.
if (samples.num_columns() != coeffs.num_rows())
return FALSE;
multiply(samples,coeffs,res);
return TRUE;
}
int stepwise_ols(const EST_FMatrix &X,
const EST_FMatrix &Y,
const EST_StrList &feat_names,
float limit,
EST_FMatrix &coeffs,
const EST_FMatrix &Xtest,
const EST_FMatrix &Ytest,
EST_IVector &included)
{
// Find the features that contribute to the correlation using a
// a greedy algorithm
EST_FMatrix coeffsl;
float best_score=0.0,bscore;
int i,best_feat;
int nf=1; // for nice printing of progress
for (i=1; i < X.num_columns(); i++)
{
if (!ols_stepwise_find_best(X,Y,included,coeffsl,
bscore,best_feat,Xtest,Ytest,
feat_names))
{
cerr << "OLS: stepwise failed" << endl;
return FALSE;
}
if ((bscore - (bscore * (limit/100))) <= best_score)
break;
else
{
best_score = bscore;
coeffs = coeffsl;
included[best_feat] = TRUE;
printf("FEATURE %d %s: %2.4f\n",
nf,
(const char *)feat_names.nth(best_feat),
best_score);
fflush(stdout);
nf++;
}
}
return TRUE;
}
static int ols_stepwise_find_best(const EST_FMatrix &X,
const EST_FMatrix &Y,
EST_IVector &included,
EST_FMatrix &coeffs,
float &bscore,
int &best_feat,
const EST_FMatrix &Xtest,
const EST_FMatrix &Ytest,
const EST_StrList &feat_names
)
{
EST_FMatrix coeffsl;
bscore = 0;
best_feat = -1;
int i;
for (i=0; i < included.length(); i++)
{
if (included.a_no_check(i) == FALSE)
{
float cor, rmse;
EST_FMatrix pred;
included.a_no_check(i) = TRUE;
if (!robust_ols(X,Y,included,coeffsl))
return FALSE; // failed for some reason
ols_apply(Xtest,coeffsl,pred);
ols_test(Ytest,pred,cor,rmse);
printf("tested %d %s %f best %f\n",
i,(const char *)feat_names.nth(i),cor,bscore);
if (fabs(cor) > bscore)
{
bscore = fabs(cor);
coeffs = coeffsl;
best_feat = i;
}
included.a_no_check(i) = FALSE;
}
}
return TRUE;
}
int ols_test(const EST_FMatrix &real,
const EST_FMatrix &predicted,
float &correlation,
float &rmse)
{
// Others probably want this function too
// return correlation and RMSE for col 0 in real and predicted
int i;
float p,r;
EST_SuffStats x,y,xx,yy,xy,se,e;
double error;
double v1,v2,v3;
if (real.num_rows() != predicted.num_rows())
return FALSE; // can't do this
for (i=0; i < real.num_rows(); i++)
{
r = real(i,0);
p = predicted(i,0);
x += p;
y += r;
error = p-r;
se += error*error;
e += fabs(error);
xx += p*p;
yy += r*r;
xy += p*r;
}
rmse = sqrt(se.mean());
v1 = xx.mean()-(x.mean()*x.mean());
v2 = yy.mean()-(y.mean()*y.mean());
v3 = v1*v2;
if (v3 <= 0)
{ // happens when there's very little variation in x
correlation = 0;
rmse = se.mean();
return FALSE;
}
// Pearson's product moment correlation coefficient
correlation = (xy.mean() - (x.mean()*y.mean()))/ sqrt(v3);
// I hate to have to do this but it is necessary.
// When the the variation of X is very small v1*v2 approaches
// 0 (the negative and equals case is caught above) but that
// may not be enough when v1 or v2 are very small but positive.
// So I catch it here. If I knew more math I'd be able to describe
// this better but the code would remain the same.
if ((correlation <= 1.0) && (correlation >= -1.0))
return TRUE;
else
{
correlation = 0;
return FALSE;
}
}
EST_FMatrix
Definition EST_FMatrix.h:56
EST_TList::nth
T & nth(int n)
return the Nth value
Definition EST_TList.h:139
EST_TMatrix::num_columns
int num_columns() const
return number of columns
Definition EST_TMatrix.h:181
EST_TMatrix::a_no_check
INLINE const T & a_no_check(int row, int col) const
const access with no bounds check, care recommend
Definition EST_TMatrix.h:184
EST_TMatrix::num_rows
int num_rows() const
return number of rows
Definition EST_TMatrix.h:179
EST_TSimpleMatrix::resize
void resize(int rows, int cols, int set=1)
resize matrix
Definition EST_TSimpleMatrix.cc:86
EST_TSimpleVector::resize
void resize(int n, int set=1)
resize vector
Definition EST_TSimpleVector.cc:70
EST_TVector::length
INLINE int length() const
number of items in vector.
Definition EST_TVector.h:252
EST_TVector::a_no_check
INLINE const T & a_no_check(int n) const
read-only const access operator: without bounds checking
Definition EST_TVector.h:257