Skip to content
Snippets Groups Projects
Commit 12c60252 authored by Paul McCarthy's avatar Paul McCarthy :mountain_bicyclist:
Browse files

MNT: Include cprob/libprob.h instead of including libprob.h directly

parent 882ef572
No related branches found
No related tags found
1 merge request!7MNT: Include cprob/libprob.h instead of including libprob.h directly
Hide whitespace changes
Inline Side-by-side
Makefile
+ 3
− 3
View file @ 12c60252
......@@ -4,14 +4,14 @@ include ${FSLCONFDIR}/default.mk
PROJNAME = miscmaths
USRINCFLAGS = -I${INC_NEWMAT} -I${INC_BOOST} -I${INC_PROB}
USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_PROB}
USRINCFLAGS = -I${INC_NEWMAT} -I${INC_BOOST} -I${INC_CPROB}
USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_CPROB}
USRCXXFLAGS = -std=c++11
OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o minimize.o cspline.o sparse_matrix.o sparsefn.o rungekutta.o nonlin.o bfmatrix.o Simplex.o SpMatMatrices.o
LIBS = -lutils -lnewmat -lprob -lm
LIBS = -lutils -lcprob -lm
# The target "all" should ALWAYS be provided
# typically it will just be another target name
......
base2z.cc
+ 10
− 40
View file @ 12c60252
......@@ -8,7 +8,7 @@
#include <cmath>
#include "base2z.h"
#include "libprob.h"
#include "cprob/libprob.h"
namespace MISCMATHS {
......@@ -19,7 +19,7 @@ namespace MISCMATHS {
return MISCMATHS::lgam(v)+MISCMATHS::lgam(w)-MISCMATHS::lgam(v+w);
}
float Base2z::logp2largez(float logp)
float Base2z::logp2largez(float logp)
{
// Large Z extrapolation routine for converting log(p) to Z values
// written by Mark Jenkinson, March 2000
......@@ -27,11 +27,11 @@ namespace MISCMATHS {
// Equations were derived by using integration by parts and give the
// following formulae:
// Z to log(p)
// log(p) = -1/2*z*z - 1/2*log(2*pi) - log(z)
// log(p) = -1/2*z*z - 1/2*log(2*pi) - log(z)
// + log(1 - 1/(z*z) + 3/(z*z*z*z))
// this equation 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)
// 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
// The equation is accurate to 1 part in 10^3 for Z>3.12 (3 iterations)
......@@ -39,25 +39,25 @@ namespace MISCMATHS {
static const float pi = 3.141592653590;
static const float log2pi = log(2*pi);
float z0, zn;
// iteratively solve for z given log p
float b = -2*logp - log2pi;
float b = -2*logp - log2pi;
z0 = sqrt(b);
zn = z0;
for (int m=1; m<=3; m++) {
// zn = sqrt(b + 2*log(1/zn - 1/(zn*zn*zn) + 3/(zn*zn*zn*zn*zn)));
zn = sqrt(b + 2*log(((3/(zn*zn) - 1)/(zn*zn) + 1)/zn) );
}
return zn;
}
float Base2z::convertlogp2z(float logp)
float Base2z::convertlogp2z(float logp)
{
// logp must be the *natural* logarithm of p, not base 10
float z = 0.0;
if(!issmalllogp(logp)) {
z = MISCMATHS::ndtri(exp(logp));
}
......@@ -69,33 +69,3 @@ namespace MISCMATHS {
}
}
f2z.cc
+ 24
− 54
View file @ 12c60252
......@@ -10,8 +10,8 @@
#include "f2z.h"
#include "utils/log.h"
#include "utils/tracer_plus.h"
#include "cprob/libprob.h"
#include <stdexcept>
#include "libprob.h"
using namespace NEWMAT;
using namespace Utilities;
......@@ -19,7 +19,7 @@ using namespace Utilities;
namespace MISCMATHS {
F2z* F2z::f2z = NULL;
float F2z::largef2logp(float f, int d1, int d2)
{
Tracer_Plus ts("F2z::largef2logp");
......@@ -35,12 +35,12 @@ namespace MISCMATHS {
cerr << "f cannot be zero or negative!" << endl;
return 0.0;
}
if (d1<=0 || d2<=0) {
if (d1<=0 || d2<=0) {
cerr << "DOFs cannot be zero or negative!" << endl;
return 0.0;
return 0.0;
}
double alpha=d1/(double)d2;
double m=(d1+d2)/2.0;
double n=(1-d1/2.0);
......@@ -51,7 +51,7 @@ namespace MISCMATHS {
double top = 1.0;
double bot = n+m-1;
double iter = 0.0;
// cerr << "logbeta(d2/2.0,d1/2.0)=" << logbeta(d2/2.0,d1/2.0) << endl;
// cerr << "loggam = " << loggam << endl;
// cerr << "n = " << n << endl;
......@@ -60,7 +60,7 @@ namespace MISCMATHS {
for(int i = 1; i <= N; i++)
{
// cerr << "i=" << i;
iter = iter + top* ( std::pow( f,float(-(n+i-1)) ) / ( std::pow(alpha,double(i))*bot ) );
iter = iter + top* ( std::pow( f,float(-(n+i-1)) ) / ( std::pow(alpha,double(i))*bot ) );
top = top*(n-1+i)*(-1);
bot = bot*(n+m-1+i);
// cerr << "iter=" << iter;
......@@ -78,15 +78,15 @@ namespace MISCMATHS {
}
bool F2z::islargef(float f, int d1, int d2, float &logp) {
if(f > 2.0 && d1>1)
{
try
{
logp=largef2logp(f,d1,d2);
logp=largef2logp(f,d1,d2);
}
catch(Exception& p_excp)
catch(Exception& p_excp)
{
cerr << "Negative iter in F2z::largef2logp" << endl;
return false;
......@@ -102,7 +102,7 @@ namespace MISCMATHS {
return (logp < -14.5);
}
float F2z::convert(float f, int d1, int d2)
float F2z::convert(float f, int d1, int d2)
{
Tracer_Plus ts("F2z::convert");
......@@ -128,18 +128,18 @@ namespace MISCMATHS {
dof2 = p_dof2;
ComputeFStats(p_fs,p_dof1,dof2,p_zs);
}
void F2z::ComputeFStats(const ColumnVector& p_fs, int p_dof1, const ColumnVector& p_dof2, ColumnVector& p_zs)
{
Tracer_Plus ts("F2z::ComputeFStats");
int numTS = p_fs.Nrows();
p_zs.ReSize(numTS);
F2z& f2z = F2z::getInstance();
for(int i = 1; i <= numTS; i++)
{
{
if (p_fs(i) > 0.0)
{
......@@ -148,25 +148,25 @@ namespace MISCMATHS {
// cerr << ",p_dof1=" << p_dof1;
// cerr << ",p_dof2=" << p_dof2(i) << endl;
p_zs(i) = f2z.convert(p_fs(i),int(p_dof1),int(p_dof2(i)));
p_zs(i) = f2z.convert(p_fs(i),int(p_dof1),int(p_dof2(i)));
}
else
{
p_zs(i) = 0.0;
}
}
}
}
void F2z::ComputeFStats(const ColumnVector& p_fs, const ColumnVector& p_dof1, const ColumnVector& p_dof2, ColumnVector& p_zs)
{
Tracer_Plus ts("F2z::ComputeFStats");
int numTS = p_fs.Nrows();
p_zs.ReSize(numTS);
F2z& f2z = F2z::getInstance();
for(int i = 1; i <= numTS; i++)
{
{
if (p_fs(i) > 0.0)
{
......@@ -175,43 +175,13 @@ namespace MISCMATHS {
// cerr << ",p_dof1=" << p_dof1;
// cerr << ",p_dof2=" << p_dof2(i) << endl;
p_zs(i) = f2z.convert(p_fs(i),int(p_dof1(i)),int(p_dof2(i)));
p_zs(i) = f2z.convert(p_fs(i),int(p_dof1(i)),int(p_dof2(i)));
}
else
{
p_zs(i) = 0.0;
}
}
}
}
}
}
miscprob.h
+ 15
− 21
View file @ 12c60252
......@@ -13,28 +13,28 @@
#define __miscprob_h
#include "miscmaths.h"
#include "libprob.h"
#include "cprob/libprob.h"
#include "stdlib.h"
using namespace NEWMAT;
namespace MISCMATHS {
// ReturnMatrix betarnd(const int dim1, const int dim2,
// const float a, const float b);
// ReturnMatrix betarnd(const int dim1, const int dim2,
// const float a, const float b);
ReturnMatrix betapdf(const RowVector& vals,
const float a, const float b);
ReturnMatrix betapdf(const RowVector& vals,
const float a, const float b);
ReturnMatrix unifrnd(const int dim1 = 1, const int dim2 = -1,
ReturnMatrix unifrnd(const int dim1 = 1, const int dim2 = -1,
const float start = 0, const float end = 1);
ReturnMatrix normrnd(const int dim1 = 1, const int dim2 = -1,
ReturnMatrix normrnd(const int dim1 = 1, const int dim2 = -1,
const float mu = 0, const float sigma = 1);
// returns nsamps*nparams matrix:
ReturnMatrix mvnrnd(const RowVector& mu, const SymmetricMatrix& covar, int nsamp = 1);
float mvnpdf(const RowVector& vals, const RowVector& mu, const SymmetricMatrix& covar);
float bvnpdf(const RowVector& vals, const RowVector& mu, const SymmetricMatrix& covar);
......@@ -44,7 +44,7 @@ namespace MISCMATHS {
ReturnMatrix normpdf(const RowVector& vals, const float mu = 0, const float var = 1);
ReturnMatrix normpdf(const RowVector& vals, const RowVector& mus,
ReturnMatrix normpdf(const RowVector& vals, const RowVector& mus,
const RowVector& vars);
ReturnMatrix normcdf(const RowVector& vals, const float mu = 0, const float var = 1);
......@@ -53,13 +53,13 @@ namespace MISCMATHS {
ReturnMatrix gammacdf(const RowVector& vals, const float mu = 0, const float var = 1);
// ReturnMatrix gammarnd(const int dim1, const int dim2,
// ReturnMatrix gammarnd(const int dim1, const int dim2,
// const float a, const float b);
// returns n! * n matrix of all possible permutations
ReturnMatrix perms(const int n);
class Mvnormrandm
{
public:
......@@ -76,14 +76,14 @@ namespace MISCMATHS {
covarw = sqrt(eig_val)*eig_vec.t();
}
ReturnMatrix next(int nsamp = 1) const
ReturnMatrix next(int nsamp = 1) const
{
Matrix ret = ones(nsamp, 1)*mu + normrnd(nsamp,mu.Ncols())*covarw;
ret.Release();
return ret;
}
ReturnMatrix next(const RowVector& pmu, int nsamp = 1)
ReturnMatrix next(const RowVector& pmu, int nsamp = 1)
{
mu=pmu;
......@@ -106,7 +106,7 @@ namespace MISCMATHS {
covarw = sqrt(eig_val)*eig_vec.t();
}
private:
private:
RowVector mu;
SymmetricMatrix covar;
......@@ -116,9 +116,3 @@ namespace MISCMATHS {
};
}
#endif
quick.cc
+ 4
− 6
View file @ 12c60252
......@@ -15,7 +15,6 @@
#include "miscmaths.h"
#include "t2z.h"
//#include "libprob.h"
using namespace MISCMATHS;
......@@ -24,21 +23,20 @@ int main(int argc, char *argv[])
try{
Matrix X = read_vest("/usr/people/woolrich/matlab/vbbabe/data/design2.mat").t();
ColumnVector Y = read_vest("/usr/people/woolrich/matlab/vbbabe/data/sdf2.mat").t();
ColumnVector m_B;
SymmetricMatrix ilambda_B;
glm_vb(X, Y, m_B, ilambda_B, 30);
write_ascii_matrix(m_B,"/usr/people/woolrich/matlab/vbbabe/data/m_B");
}
catch(Exception p_excp)
catch(Exception p_excp)
{
cerr << p_excp.what() << endl;
}
return 0;
}
t2z.cc
+ 20
− 50
View file @ 12c60252
......@@ -10,7 +10,7 @@
#include "t2z.h"
#include "newmat.h"
#include "utils/tracer_plus.h"
#include "libprob.h"
#include "cprob/libprob.h"
using namespace NEWMAT;
using namespace Utilities;
......@@ -19,8 +19,8 @@ namespace MISCMATHS {
T2z* T2z::t2z = NULL;
Z2t* Z2t::z2t = NULL;
float Z2t::convert(float z, int dof)
float Z2t::convert(float z, int dof)
{
float t = 0.0;
......@@ -37,9 +37,9 @@ namespace MISCMATHS {
float T2z::larget2logp(float t, int dof)
{
// static float logbeta[] = { 1.144729885849, 0.693147180560,
// 0.451582705289, 0.287682072452,
// 0.451582705289, 0.287682072452,
// 0.163900632838, 0.064538521138,
// -0.018420923956, -0.089612158690,
// -0.018420923956, -0.089612158690,
// -0.151952316581, -0.207395194346 } ;
//static const float pi = 3.141592653590;
......@@ -54,13 +54,13 @@ namespace MISCMATHS {
// 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))
// + 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(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)
// 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)
......@@ -70,12 +70,12 @@ namespace MISCMATHS {
if (t<0) {
return larget2logp(-t,dof);
}
float logp, lbeta;
if (dof<=0) {
if (dof<=0) {
cerr << "DOF cannot be zero or negative!" << endl;
return 0.0;
return 0.0;
}
float n = (float) dof;
......@@ -87,7 +87,7 @@ namespace MISCMATHS {
//} 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;
......@@ -98,7 +98,7 @@ namespace MISCMATHS {
}
bool T2z::islarget(float t, int dof, float &logp) {
// aymptotic formalae are valid if
// 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
......@@ -118,7 +118,7 @@ namespace MISCMATHS {
float z = 0.0, logp=0.0;
double p(0);
if(!islarget(t,dof,logp)) {
// cerr << "t = " << t << endl;
p = MISCMATHS::stdtr(dof, t);
......@@ -127,7 +127,7 @@ namespace MISCMATHS {
z = MISCMATHS::ndtri(p);
}
else {
z = logp2largez(logp);
// cerr<<endl<<"logp="<<logp<<endl;
......@@ -143,16 +143,16 @@ namespace MISCMATHS {
}
float T2z::converttologp(float t, int dof)
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
// but t to logp is not available for negative t
// so just hardcode it to be -1e-12
logp=-1e-12;
}
......@@ -187,7 +187,7 @@ namespace MISCMATHS {
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;
}
......@@ -228,7 +228,7 @@ namespace MISCMATHS {
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;
}
......@@ -238,33 +238,3 @@ namespace MISCMATHS {
}
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment