From 3719572c81ab3ee4c5ea87ddaff1db049fe0854e Mon Sep 17 00:00:00 2001
From: Mark Jenkinson <mark@fmrib.ox.ac.uk>
Date: Tue, 21 Aug 2001 18:52:05 +0000
Subject: [PATCH] Allowing dimensionality to change - also fixed a small maths
error
---
infer.cc | 70 +++++++++++++++++++++++++++++++++++++++++++-------------
1 file changed, 54 insertions(+), 16 deletions(-)
diff --git a/infer.cc b/infer.cc
index cff380c..bf4ad24 100644
--- a/infer.cc
+++ b/infer.cc
@@ -28,23 +28,28 @@ Infer::Infer(float udLh, float ut, unsigned int uV) {
t = ut;
V = uV;
if (V<=0.0) V=1.0;
+ // dimensionality
+ D = 3.0;
+ // if (zsize <= 1) D = 2.0; // to be done by calling program
+ // NB: the (sqr(t) -1) is previous D=3 version (from where??)
if (fabs(t)<13.0) {
- Em_ = V * pow(2*M_PI,-2) * dLh * (sqr(t) - 1) * exp(-sqr(t)/2.0);
+ Em_ = V * pow(2*M_PI,-(D+1)/2) * dLh * pow((sqr(t) - 1), (D-1)/2) *
+ exp(-sqr(t)/2.0);
} else {
Em_ = 0.0; // underflowed exp()
}
if (fabs(t)<8.0) {
- B_ = pow((gamma(2.5)*Em_)/(V*(0.5 + 0.5*erf(-t/sqrt(2)))),(2.0/3.0));
+ B_ = pow((gamma(1.0+D/2.0)*Em_)/(V*(0.5 + 0.5*erf(-t/sqrt(2)))),(2.0/D));
} else {
- // the large t approximation
- float a1 = V * dLh / (4 * M_PI * M_PI);
- float a3 = pow((gamma(2.5) / V ),(2.0/3.0));
+ // the large t approximation (see appendix below)
+ float a1 = V * dLh * pow(2*M_PI,-(D+1)/2);
+ float a3 = pow((gamma(1+D/2.0) / V ),(2.0/D));
float tsq = t*t;
- float Em_q = a1 * (2.0 * M_PI) * (tsq - 1.0)*t
- / ( 1.0 - 1.0/tsq + 3.0/(tsq*tsq)) ;
- B_ = a3 * pow(Em_q,(2.0/3.0));
+ float c = pow(2*M_PI,-1.0/2.0) * t / ( 1.0 - 1.0/tsq + 3.0/(tsq*tsq)) ;
+ float Em_q = a1 * pow((tsq - 1.0),(D-1)/2) * c;
+ B_ = a3 * pow(Em_q,(2.0/D));
}
@@ -54,10 +59,10 @@ Infer::Infer(float udLh, float ut, unsigned int uV) {
float Infer::operator() (unsigned int k) {
// ideally returns the following:
- // return 1 - exp(-Em_ * exp(-B_ * pow( k , 2.0 / 3.0)));
+ // return 1 - exp(-Em_ * exp(-B_ * pow( k , 2.0 / D)));
// but in practice must be careful about ranges
// Assumes that exp(+/-87) => 1e+/38 is OK for floats
- float arg1 = -B_ * pow(k , 2.0 / 3.0);
+ float arg1 = -B_ * pow(k , 2.0 / D);
if (fabs(arg1)<87.0) {
float exp1 = exp(arg1);
float arg2 = -Em_ * exp1;
@@ -66,16 +71,49 @@ float Infer::operator() (unsigned int k) {
return 1.0 - exp2;
}
}
- // now must invoke the approximation instead
- float a1 = V * dLh / (4 * M_PI * M_PI);
- float a3 = pow((gamma(2.5) / V ),(2.0/3.0));
+ // now must invoke the approximation instead (see appendix below)
+ float a1 = V * dLh * pow(2*M_PI,-(D+1)/2);
+ float a3 = pow((gamma(1.0+D/2.0) / V ),(2.0/D));
float tsq = t*t;
- float c = 1.0 / ( 1.0 - 1.0/tsq + 3.0/(tsq*tsq));
- float arg = -tsq*(0.5 + a3*pow(a1,2.0/3.0)*pow(2.0*M_PI,1.0/3.0)*c*pow(k,2.0/3.0));
+ float c = pow(2*M_PI,-1.0/2.0) * t / ( 1.0 - 1.0/tsq + 3.0/(tsq*tsq));
+ float Em_q = a1 * pow((tsq - 1.0),(D-1)/2) * c;
+ float beta = a3 * pow(Em_q , 2.0/D);
+ float arg = -tsq*0.5 - beta*pow(k,2.0/D));
float exp1 = 0.0;
if (arg>-87.0) {
exp1 = exp(arg);
}
- float p = a1*(tsq-1.0)*exp1;
+ float p = a1*pow((tsq-1.0),(D-1)/2)*exp1;
return p;
}
+
+
+
+// MATHEMATICAL APPENDIX
+
+/*
+
+The formulas that need to be calculated are:
+(1) E_m = V * dLh * (2*pi)^(-(D+1)/2) * (t^2 -1)^((D-1)/2) * exp(-t^2 /2)
+(2) Beta = (Gamma(D/2+1)/V * E_m / Phi(-t) )^(2/D)
+(3) p = 1 - exp( - E_m * exp(-Beta*k^(2/D)))
+
+where Phi(-t) = Gaussian cumulant = (1/2 + 1/2*erf(-t/sqrt(2)))
+
+These are approximated by:
+
+(2a) Beta = (Gamma(D/2+1)/V)^(2/D) * (Em1)^(2/D) * Ct^(2/D)
+where Em1 = V * dLh * (2*pi)^(-(D+1)/2) * (t^2 -1)^((D-1)/2)
+ Ct = (2*pi)^(-1/2) * t / ( 1.0 - 1.0/t^2 + 3.0/t^4 )
+ which approximates ( exp(-t^2 /2) / Phi(-t) )^(2/D)
+ using 1/2 - 1/2*erf(t/sqrt(2)) = (2*pi)^(1/2) * exp(-t^2 /2) * ...
+ (1-1/t^2+3/t^4) / t
+ (this is derived in TR00MJ1)
+
+and
+
+(3a) p = Em1 * exp( -t^2 /2 - beta * k^(2/D))
+using the approximation (2a) for beta
+
+ */
+
--
GitLab