-
Mark Jenkinson authoredMark Jenkinson authored
t2z.cc 6.12 KiB
/* t2z.cc
Mark Woolrich & Mark Jenkinson, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
#include <cmath>
#include "t2z.h"
#include "newmat.h"
#include "utils/tracer_plus.h"
#include "libprob.h"
using namespace NEWMAT;
using namespace Utilities;
namespace MISCMATHS {
T2z* T2z::t2z = NULL;
Z2t* Z2t::z2t = NULL;
float Z2t::convert(float z, int dof)
{
float t = 0.0;
if(z>8)
throw Exception("z is too large to convert to t");
double p = MISCMATHS::ndtr(z);
cerr << "p = " << p << endl;
t = MISCMATHS::stdtri(dof,p);
return t;
}
float T2z::larget2logp(float t, int dof)
{
// static float logbeta[] = { 1.144729885849, 0.693147180560,
// 0.451582705289, 0.287682072452,
// 0.163900632838, 0.064538521138,
// -0.018420923956, -0.089612158690,
// -0.151952316581, -0.207395194346 } ;
//static const float pi = 3.141592653590;
// static const float log2pi = log(2*pi);
// Large T extrapolation routine for converting T to Z values
// written by Mark Jenkinson, March 2000
//
// It does T to Z via log(p) rather than p, since p becomes very
// small and underflows the arithmetic
// Equations were derived by using integration by parts and give the
// following formulae:
// (1) T to log(p) NB: n = DOF
// log(p) = -1/2*log(n) - log(beta(n/2,1/2)) - (n-1)/2*log(1+t*t/n)
// + log(1 - (n/(n+2))/(t*t) + 3*n*n/((n+2)*(n+4)*t*t*t*t))
// (2) Z to log(p)
// log(p) = -1/2*z*z - 1/2*log(2*pi) - log(z)
// + log(1 - 1/(z*z) + 3/(z*z*z*z))
// equation (2) is then solved by the recursion:
// z_0 = sqrt(2*(-log(p) - 1/2*log(2*pi)))
// z_{n+1} = sqrt(2*(-log(p) - 1/2*log(2*pi) - log(z_n)
// + log(1 - 1/(zn*zn) + 3/(zn*zn*zn*zn))
// In practice this recursion is quite accurate in 3 to 5 iterations
// Equation (1) is accurate to 1 part in 10^3 for T>7.5 (any n)
// Equation (2) is accurate to 1 part in 10^3 for Z>3.12 (3 iterations)
if (t<0) { return larget2logp(-t,dof);
}
float logp, lbeta;
if (dof<=0) {
cerr << "DOF cannot be zero or negative!" << endl;
return 0.0;
}
float n = (float) dof;
// complete Beta function
lbeta = this->logbeta(1/2.0,n/2.0);
//if (dof<=10) {
//lbeta = logbeta[dof-1];
//} else {
//lbeta = log2pi/2 - log(n)/2 + 1/(4*n);
//}
// log p from t value
// logp = log( (1 - n/((n+2)*t*t) + 3*n*n/((n+2)*(n+4)*t*t*t*t))/(sqrt(n)*t))
// - ((n-1)/2)*log(1 + t*t/n) - lbeta;
logp = log(( (3*n*n/((n+2)*(n+4)*t*t) - n/(n+2))/(t*t) + 1)/(sqrt(n)*t))
- ((n-1)/2)*log(1 + t*t/n) - lbeta;
return logp;
}
bool T2z::islarget(float t, int dof, float &logp) {
// aymptotic formalae are valid if
// log(p) < -14.5 (derived from Z-statistic approximation error)
// For dof>=15, can guarantee that log(p)>-33 (p > 1e-14) if T<7.5
// and so in this region use conventional means, not asymptotic
if ((dof>=15) && (fabs(t)<7.5)) { return false; }
logp=larget2logp(t,dof);
if (dof>=15) return true; // force asymptotic calc for all T>=7.5, D>=15
return issmalllogp(logp);
}
bool T2z::issmalllogp(float logp) {
// aymptotic formula accurate to 1 in 10^3 for Z>4.9
// which corresponds to log(p)=-14.5
return (logp < -14.5);
}
float T2z::convert(float t, int dof) {
float z = 0.0, logp=0.0;
if(!islarget(t,dof,logp)) {
// cerr << "t = " << t << endl;
double p = MISCMATHS::stdtr(dof, t);
//cerr << "p = " << p << endl;
z = MISCMATHS::ndtri(p);
}
else {
z = logp2largez(logp);
// cerr<<endl<<"logp="<<logp<<endl;
if (t<0) z=-z;
}
return z;
}
float T2z::converttologp(float t, int dof)
{
float logp=0.0;
if(!islarget(t,dof,logp)) {
logp = log(1-MISCMATHS::stdtr(dof, t));
}
else if(t<0) {
// t < 0 and abs(t) is large enough to require asymptotic approx.
// but t to logp is not available for negative t
// so just hardcode it to be -1e-12
logp=-1e-12;
}
// cerr << "logp = " << logp << endl;
// cerr << "exp(logp) = " << std::exp(logp) << endl;
return logp;
}
void T2z::ComputePs(const ColumnVector& p_vars, const ColumnVector& p_cbs, int p_dof, ColumnVector& p_ps)
{
Tracer ts("T2z::ComputePs");
int numTS = p_vars.Nrows();
T2z& t2z = T2z::getInstance();
p_ps.ReSize(numTS);
for(int i = 1; i <= numTS; i++)
{
//cerr << "var = " << p_vars(i) << " at index "<< i << endl;
//cerr << "cb = " << p_cbs(i) << " at index "<< i << endl;
if ( (p_vars(i) != 0.0) && (p_cbs(i) != 0.0) )
{
if(p_vars(i) < 0.0)
{
//cerr << "var = " << p_vars(i) << " at index "<< i << endl;
p_ps(i) = 0.0;
}
else
{
p_ps(i) = t2z.converttologp(p_cbs(i)/sqrt(p_vars(i)),p_dof);
//if(p_zs(i) == 0.0)
//cerr << " at index " << i << endl;
}
}
else
p_ps(i) = 0.0;
}
}
void T2z::ComputeZStats(const ColumnVector& p_vars, const ColumnVector& p_cbs, int p_dof, ColumnVector& p_zs)
{
ColumnVector dof = p_vars;
dof = p_dof;
ComputeZStats(p_vars,p_cbs,dof,p_zs);
}
void T2z::ComputeZStats(const ColumnVector& p_vars, const ColumnVector& p_cbs, const ColumnVector& p_dof, ColumnVector& p_zs)
{
Tracer ts("T2z::ComputeStats");
int numTS = p_vars.Nrows();
T2z& t2z = T2z::getInstance();
p_zs.ReSize(numTS);
for(int i = 1; i <= numTS; i++)
{
//cerr << "var = " << p_vars(i) << " at index "<< i << endl;
//cerr << "cb = " << p_cbs(i) << " at index "<< i << endl;
if ( (p_vars(i) != 0.0) && (p_cbs(i) != 0.0) )
{
if(p_vars(i) < 0.0)
{
//cerr << "var = " << p_vars(i) << " at index "<< i << endl;
p_zs(i) = 0.0;
}
else
{
p_zs(i) = t2z.convert(p_cbs(i)/sqrt(p_vars(i)),int(p_dof(i)));
//if(p_zs(i) == 0.0)
//cerr << " at index " << i << endl;
}
}
else
p_zs(i) = 0.0;
}
}
}