Colin Champion, 17 Oct 2022
This page provides a short C function to integrate a bivariate Gaussian over the region defined by at most 3 linear inequalities (or over a disc).
gaussint : software : test : download : index
Suppose that we have a standard bivariate Gaussian, and that we wish to integrate the Gaussian over the region satisfying three inequalities. This integral can only be computed numerically (and approximately). Its value is given to good precision by the call
result = gaussint(lt,rhs,nlt,tol,maxeval) ;
where:
The return value is an xy structure (as defined in memory.h) whose x field is the desired integral and whose y field is a condition number, which will be <1 (and ≥0) if the integral converged, and >1 otherwise.
To be specific, gaussint attempts to reduce the estimated error to below tol (relatively or absolutely), and the condition number is the ratio of the estimated error to tol on completion. Thus something slightly greater than 1 doesn’t matter, but when convergence is unattainable the condition number is likely to be huge (and the returned integral wholly unreliable).
A similar functions gaussdisc is provided to integrate a standard Gaussian over a disc.
If inequality 0 is “2x – 3y > 4” then lt[0][0] will be –2, lt[0][1] will be 3, and rhs[0] will be –4 (the signs of all quantities having been flipped to convert a “>” into a “<”).
All floating point numbers are doubles.
gaussint may apply a relative or an absolute tolerance. If the value you supply is positive, eg. 10–6, then gaussint will try to find an answer to within one part in a million. If you supply a negative tolerance, eg. –10–6, then gaussint will try to find an answer which is correct to within 0.000001. In the tails of the distribution, all results will be 0 to within any reasonable absolute tolerance.
If you try to push the tolerance too far, you’ll be let down by the limited precision of the Gaussian tail areas computed within gaussint (which I think use coefficients from “Numerical Recipes in C”). If you can find better coefficients, you may be able to get more accurate results.
result = gaussint(lt,rhs,nlt,tol,maxeval) ;
as above. There is also a simplified call
integral = gaussint(lt,rhs,nlt) ;
in which the integral is performed with an absolute tolerance of 10–7 and at most 20000 function evaluations. The simple call will crash out if the tolerance isn’t attained, so you can trust any result returned.
A separate function is provided for the relatively simple case of integrating over a triangle:
result = gausstri(A,B,C,tol,maxeval) ;
where A, B and C are xy’s containing the coordinates of the vertices. This too has a simplified call:
integral = gausstri(A,B,C) ;
as before.
An additional function integrates a Gaussian over a disc and outside the disc.
result = gaussdisc(d,r,tol,maxeval) ; integral = gaussdisc(d,r) ; result = gaussexdisc(d,r,tol,maxeval) ; integral = gaussexdisc(d,r) ;
r is the radius of the disc and d is the distance of its centre from the origin. gaussdisc integrates the interior and gaussexdisc the exterior. Obviously the sum or the two integrals is 1, but if you try to obtain one integral from the other, and the other is close to 1, then you will lose precision (both because of roundoff error, and because a relative tolerance for integration is relative to the wrong value). A rule of thumb is to use gaussdisc when r<d+1/(2d+1) and to use gaussexdisc when r>d+1/(2d+1).
If you want to integrate round the circumference, write
double gausscirc(double d,double r)
{ return r * exp(-(d*d+r*r)/2) * besseli(0,r*d) ; }
where the modified Bessel function Iν(x) is provided along with my spatial solver.
The limiting cases of these functions are occasionally of interest.
| ● gaussdisc(0,r) = 1 – exp(–r2/2) | ● gausscirc(0,r) = r exp(–r2/2) | ||
| ● gaussdisc(d,ε) = ε2 exp(–d2/2) / 2 | ● gausscirc(d,ε) = ε exp(–d2/2) | ||
| ● gaussdisc(∞,r) = 0 | ● gausscirc(D,r) = √(r/(2πD)) exp(–(D–r)2/2) | as D → ∞ | |
| ● gaussdisc(d,R) = 1 – √(R/(2πd)) exp(–(R–d)2/2) / (R–d) [d>0] | ● gausscirc(d,R) = √(R/(2πd)) exp(–(R–d)2/2) [d>0] | as R → ∞ | |
| ● gaussdisc(0,∞) = undefined | ● gausscirc(0,∞) = undefined | ||
| ● gaussdisc(∞,∞) = undefined | ● gausscirc(∞,∞) = undefined. |
• swap • radians • euclid2 • euclid • gausslog1 • gauss1 • gausslog2 • gauss2 • rotate • solv2 • gausstri • gaussint • discint • gaussdisc • gaussexdisc
#include "memory.h" #include <math.h> #include "quadrate.h" static void swap(double *&a,double *&b) { double *x=a ; a = b ; b = x ; } static double radians() { return 2 * 3.1415926535897932384626433832795029 ; } xy linmax(double (*f)(double),xy A,xy B,xy C,double tol) ; static double radians(double x,double qlo) { static double pi=radians()/2 ; while(x>=qlo+2*pi) x -= 2*pi*(int)((x-qlo)/(2*pi)) ; while(x<qlo) x += 2*pi*(int)((qlo+2*pi-x)/(2*pi)) ; return x ; } static double radians(double x) { return radians(x,0) ; } static double euclid2(xy A,xy B) { return (A.x-B.x)*(A.x-B.x) + (A.y-B.y)*(A.y-B.y) ; } static double euclid(xy A,xy B) { return sqrt(euclid2(A,B)) ; } static double euclid2(xy A) { return euclid2(A,xy()) ; } static double euclid(xy A) { return euclid(A,xy()) ; } /* ----------------------------- objective functions ------------------------ */ double rtlnorm(double),ltlnorm(double),logrtlnorm(double),logltlnorm(double) ; double log1minus(double) ; static double pi=radians()/2,qpi=1/sqrt(2*pi),qscl,logqscl ; struct piecewise { double m0,m1,c0,c1,x0,x1 ; } ; static piecewise lo,hi ; static xy buf[66] ; static int nbuf ; static double gausslog1(double u) { double x,y ; if(u>=1) throw up ("%.1e passed to gausslog1",u) ; x = lo.x0 + qscl*atanh(u) ; if(x>=lo.x1) y = lo.m1*x + lo.c1 ; else y = lo.m0*x + lo.c0 ; y = logqscl + logrtlnorm(y) - x*x/2 - log1minus(u*u) ; buf[nbuf%64] = xy(u,y) ; nbuf += 1 ; return y ; } static void gauss1(double u,double *v) { if(u>=1) v[0] = 0 ; else v[0] = exp(gausslog1(u)) ; } static double gausslog2(double x) { double ylo=lo.c0+lo.m0*x,yhi=hi.c0+hi.m0*x,u,v ; if(ylo+yhi>0) { u = logrtlnorm(ylo) ; v = logrtlnorm(yhi) ; } else { u = logltlnorm(yhi) ; v = logltlnorm(ylo) ; } u = u - x*x/2 + log1minus(exp(v-u)) ; buf[nbuf%64] = xy(x,u) ; nbuf += 1 ; return u ; } static void gauss2(double x,double *y) { y[0] = exp(gausslog2(x)) ; } static xy rotate(xy X,double theta) { double r=euclid(X),phi=atan2(X.y,X.x) ; return xy(r*cos(theta+phi),r*sin(theta+phi)) ; } static xy solv2(double *le0,double r0,double *le1,double r1) { double q = 1 / (le0[0]*le1[1]-le0[1]*le1[0]) ; return xy( q * ( le1[1]*r0 - le0[1]*r1 ) , q * ( le0[0]*r1 - le1[0]*r0 ) ) ; } /* ---------------------------------- gausstri ------------------------------ */ xy gausstri(xy A,xy B,xy C,double tol,int maxeval) { int i ; double q,qa,lab=euclid2(A,B),lbc=euclid2(B,C),lca=euclid2(C,A),theta,**pt ; quadratereport report ; xy X,Y,Z ; // rearrange the triangle so that AB is the longest side, this side is // vertical, and B.y>A.y (so C.x may be to the left or right of A.x) if(lbc>lab&&lbc>lca) swap(A,C) ; else if(lca>lab&&lca>lbc) swap(B,C) ; theta = pi/2 - atan2(B.y-A.y,B.x-A.x) ; A = rotate(A,theta) ; B = rotate(B,theta) ; C = rotate(C,theta) ; if(A.y>B.y) swap(A,B) ; lo.x0 = A.x ; lo.x1 = C.x ; lo.m0 = (C.y-A.y) / (C.x-A.x) ; hi.m0 = (C.y-B.y) / (C.x-B.x) ; if(C.x<A.x) swap(lo.x0,lo.x1) ; lo.c0 = C.y - lo.m0*C.x ; hi.c0 = C.y - hi.m0*C.x ; // find points (X,y) such that the integrand is less at X than at Y and that // X is closer to C than is Y qa = gausslog2(A.x) ; Z = xy(A.x,qa) ; qa = exp(qa) ; nbuf = 0 ; Y = xy((A.x+C.x)/2,gausslog2) ; for(X=xy((Y.x+C.x)/2,gausslog2);X.y>Y.y;Z=Y,Y=X,X=xy((X.x+C.x)/2,gausslog2)) ; // it would be possible (and perhaps desirable), if the max is closer to Z // than to Y, to try to pin it down more closely if(Y.y>X.y&&Y.y>Z.y) Y = linmax(gausslog2,X,Y,Z,1e-4) ; q = exp(Y.y) ; if(q==0&&qa==0) return xy() ; nbuf = min(nbuf,64) ; for(i=0;i<nbuf;i++) buf[i].y = exp(buf[i].y) ; buf[nbuf] = xy(A.x,qa) ; buf[nbuf+1] = xy(C.x,0.0) ; xysort(buf,nbuf+2) ; pt = matrix(nbuf+2,2) ; for(i=0;i<nbuf+2;i++) { pt[i][0] = buf[i].x ; pt[i][1] = buf[i].y ; } report = quadraten(gauss2,1,&q,pt,nbuf+2,tol,maxeval) ; freematrix(pt) ; return xy(qpi*q,report.condition) ; } double gausstri(xy A,xy B,xy C) { A = gausstri(A,B,C,-1e-7,20000) ; if(A.y>1) throw up("convergence failure in gausstri") ; else return A.x ; } /* ---------------------------------- gaussint ------------------------------ */ xy gaussint(double **lt,double *rhs,int nlt,double tol,int maxeval) { double q,theta0,theta2,*lt0,*lt1,*lt2,r0,r1,r2,phi,eps=1e-6,**pt ; quadratereport report ; xy X,Y,Z ; int satx,saty,satz,i ; if(nlt>=1) { lt0 = lt[0] ; r0 = rhs[0] ; } if(nlt>=2) { lt1 = lt[1] ; r1 = rhs[1] ; } if(nlt>=3) { lt2 = lt[2] ; r2 = rhs[2] ; } lo.m0 = lo.m1 = lo.c0 = lo.c1 = lo.x0 = lo.x1 = qscl = logqscl = 0 ; hi = lo ; if(nlt==3) // find which of 4 poss cases: inconsistent constraints, open // segment with 2 sides, open jaw with 3 sides, triangle { X = solv2(lt1,r1,lt2,r2) ; Y = solv2(lt0,r0,lt2,r2) ; Z = solv2(lt0,r0,lt1,r1) ; satx = ( lt0[0]*X.x+lt0[1]*X.y<r0 ) ; saty = ( lt1[0]*Y.x+lt1[1]*Y.y<r1 ) ; satz = ( lt2[0]*Z.x+lt2[1]*Z.y<r2 ) ; if(satx==0&&saty==0&&satz==0) return xy() ; else if(satx&&saty&&satz) return gausstri(X,Y,Z,tol,maxeval) ; else if(satx==0&&saty==0&&satz) nlt = 2 ; else if(satx==0&&saty&&satz==0) { nlt = 2 ; lt1 = lt2 ; r1 = r2 ; } else if(satx&&saty==0&&satz==0) { nlt = 2 ; lt0 = lt2 ; r0 = r2 ; } } if(nlt==1) { q = sqrt(r0*r0 / ( lt0[0]*lt0[0] + lt0[1]*lt0[1] )) ; if(r0<0) q = rtlnorm(q) ; else q = ltlnorm(q) ; return xy(q,0.0) ; } else if(nlt==2||nlt==3) { // rotate the open jaw so that the most anticlockwise edge is vertical if(nlt==2) { lt2 = lt1 ; r2 = r1 ; Y = solv2(lt0,r0,lt2,r2) ; } else if(satx==0) { swap(X,Y) ; swap(lt0,lt1) ; swap(r0,r1) ; } else if(satz==0) { swap(Z,Y) ; swap(lt2,lt1) ; swap(r2,r1) ; } theta0 = atan2(-lt0[0],lt0[1]) ; theta2 = atan2(-lt2[0],lt2[1]) ; if(lt2[0]*cos(theta0)+lt2[1]*sin(theta0)>=0) theta0 += pi ; if(lt0[0]*cos(theta2)+lt0[1]*sin(theta2)>=0) theta2 += pi ; if(radians(theta2-theta0)>pi) { swap(theta0,theta2) ; if(nlt==3) { swap(X,Z) ; swap(lt0,lt2) ; swap(r0,r2) ; } } phi = pi/2 - theta2 ; lo.m1 = lo.m0 = tan(theta0+phi) ; Y = rotate(Y,phi) ; lo.c1 = lo.c0 = Y.y - lo.m0*Y.x ; lo.x1 = lo.x0 = Y.x ; if(nlt==3) { X = rotate(X,phi) ; Z = rotate(Z,phi) ; lo.x1 = Z.x ; lo.m0 = (Z.y-X.y) / (Z.x-X.x) ; lo.c0 = X.y - lo.m0*X.x ; } // adjust the scale factor in the change of variables so that the function // max is to the left of 0.9 - otherwise all the action may be squeezed into // a tiny sliver qscl = 1 ; logqscl = 0 ; Y = xy(0.8,gausslog1) ; Z = xy(0.9,gausslog1) ; while(Y.y<Z.y) { logqscl = log(qscl*=2) ; Y.y = gausslog1(Y.x) ; Z.y = gausslog1(Z.x) ; } // now hillclimb - first make sure we have the max bracketed X = xy(0,gausslog1) ; nbuf = 2 ; buf[0] = Y ; buf[1] = Z ; while(X.y>Y.y&&Y.x>1e-4) { Z = Y ; Y = xy(Y.x/5,gausslog1) ; } if(Y.y>X.y&&Y.y>Z.y) Y = linmax(gausslog1,X,Y,Z,1e-4) ; if(exp(Y.y)==0) return xy() ; nbuf = min(nbuf,64) ; xysort(buf,nbuf) ; pt = matrix(nbuf+2,2) ; pt[0][0] = 0 ; pt[0][1] = exp(X.y) ; for(i=0;i<nbuf;i++) { pt[1+i][0] = buf[i].x ; pt[1+i][1] = exp(buf[i].y) ; } pt[nbuf+1][0] = 1 ; pt[nbuf+1][1] = 0 ; report = quadraten(gauss1,1,&q,pt,nbuf+2,tol,maxeval) ; freematrix(pt) ; return xy(qpi*q,report.condition) ; } else if(nlt>3) throw up("gaussint accepts at most 3 inequalities") ; return xy(1,0.0) ; // nlt <= 0 } double gaussint(double **lt,double *rhs,int nlt) { xy X = gaussint(lt,rhs,nlt,-1e-7,20000) ; if(X.y>1) throw up("convergence failure in gaussint") ; else return X.x ; } /* ----------------------------------- gaussdisc ---------------------------- */ static double d,r ; static int exflag ; static void discint(double x,double *y) { double q ; if(x<0||x>=r) { y[0] = 0 ; return ; } q = sqrt(r*r-x*x) ; if(exflag==0) y[0] = qpi * exp(-x*x/2) * ( rtlnorm(d-q) - rtlnorm(d+q) ) ; else y[0] = qpi * exp(-x*x/2) * ( ltlnorm(d-q) + rtlnorm(d+q) ) ; } xy gaussdisc(double dd,double rr,double tol,int maxeval,int xf) { if(rr<=0) return xy(0.0,0.0) ; double q , **pt=matrix(3,2) ; d = dd ; r = rr ; exflag = xf ; pt[0][0] = 0 ; discint(0,&pt[0][1]) ; pt[1][0] = r/10 ; discint(r/10,&pt[1][1]) ; pt[2][0] = r ; pt[2][1] = 0 ; quadratereport report = quadraten(discint,1,&q,pt,3,tol,maxeval) ; freematrix(pt) ; if(xf) q += rtlnorm(r) ; return xy(2*q,report.condition) ; } xy gaussdisc(double dd,double rr,double tol,int maxeval) { return gaussdisc(dd,rr,tol,maxeval,0) ; } double gaussdisc(double dd,double rr) { xy X = gaussdisc(dd,rr,-1e-7,20000,0) ; if(X.y>1) throw up("convergence failure in gaussdisc") ; else return X.x ; } xy gaussexdisc(double dd,double rr,double tol,int maxeval) { return gaussdisc(dd,rr,tol,maxeval,1) ; } double gaussexdisc(double dd,double rr) { xy X = gaussdisc(dd,rr,-1e-7,20000,1) ; if(X.y>1) throw up("convergence failure in gaussdisc") ; else return X.x ; }
#include "memory.h" #include <math.h> #include "random.h" xy gaussint(double **lt,double *rhs,int nlt,double tol,int maxeval) ; // the proby of getting <n successes from N trials static double ltlbin(int n,int N,double p) { int i ; double u,v,q=1-p,r,s ; for(p/=q,s=r=1,i=1;i<n;i++) s += ( r *= p*(N+1-i)/i ) ; return s*pow(q,N) ; } int main() { int testno,ntests=100000,count,i,n=100000,sno,ex,nlt ; xy A,B,C,X ; double q,qcond,qint,**lt=matrix(3,2),*rhs=vector(3) ; for(testno=0;testno<ntests;testno++) // if(testno+1==274368) { ranset(testno) ; if((testno+1)%10000==0) printf("%d...\n",testno+1) ; i = rani(10) ; if(i<2) nlt = i ; else if(i<4) nlt = 2 ; else nlt = 3 ; for(i=0;i<nlt;i++) { lt[i][0] = gaussv() ; lt[i][1] = gaussv() ; rhs[i] = gaussv() ; } X = gaussint(lt,rhs,nlt,1e-6,20000) ; qint = X.x ; qcond = X.y ; for(count=sno=0;sno<n;sno++) { X = xy(gaussv(),gaussv()) ; for(i=0;i<nlt&<[i][0]*X.x+lt[i][1]*X.y<rhs[i];i++) ; if(i==nlt) count += 1 ; } if(qint==0||qint==1) { ex = qint?count:0 ; if(count!=ex||qcond>1) printf("%d (cond=%.2f): obs/exp = %d/%d\n",testno,X.y,count,ex) ; } else { q = (count-n*qint)/sqrt(n*qint*(1-qint)) ; if(fabs(q)>4||qcond>1) { if(qint>0&&n*qint<100&&count>n*qint) { q = 1 - ltlbin(count,n,qint) ; if(q<0.00003||qcond>1) printf("%d (cond=%.2f): proby = %.1e%s, obs/exp=%d/%.3f\n", testno,qcond,q,q<3e-7?" ***":"",count,n*qint) ; } else printf("%d (cond=%.2f): sigmage =%4.1f%s, obs/exp=%d/%d\n", testno,qcond,q,fabs(q)>5?" ***":"",count,(int)(0.5+n*qint)) ; } } } }
This generates 100 000 random sets of inequalities, calling gaussint with a fairly tight relative tolerance for each one. It generates 100 000 random samples and compares the observed number satisfying the inequalities with the expected number from the numerical integral. If there is a discrepency, or the condition number is >1, a line is printed out. I get the following results:
674 (cond=0.84): sigmage = 4.1, obs/exp=7328/6995 5190 (cond=0.87): sigmage = 4.1, obs/exp=3086/2869 6432 (cond=3869.71): sigmage =-0.0, obs/exp=0/0 9560 (cond=0.91): sigmage = 4.1, obs/exp=7392/7062 10000... 13536 (cond=0.94): sigmage =-4.1, obs/exp=9986/10378 20000... 30000... 38088 (cond=1716.24): sigmage =-0.0, obs/exp=0/0 40000... 50000... 55262 (cond=8343375000.00): sigmage =-0.0, obs/exp=0/0 60000... 70000... 79886 (cond=0.82): sigmage = 4.4, obs/exp=1725/1552 80000... 81422 (cond=4121.88): sigmage =-0.0, obs/exp=0/0 86919 (cond=0.00): sigmage =-4.1, obs/exp=99998/100000 87397 (cond=3384.28): sigmage =-0.0, obs/exp=0/0 90000... 92641 (cond=0.83): sigmage = 4.3, obs/exp=22655/22091 95144 (cond=67.07): sigmage =-0.0, obs/exp=0/0 96473 (cond=43121621.62): sigmage =-0.0, obs/exp=0/0 99182 (cond=13318600000.00): sigmage =-0.0, obs/exp=0/0
Two sorts of problem will be noticed. Firstly, occasionally there is a sigmage greater than 4 even though the condition number is <1. This is to be expected, since 4σ isn’t that rare (and the equivalent count is even less rare, given that the true distribution is binomial rather than normal. I take this into account when the expected number of occurrences is small.)
The other problem is large condition numbers. These only arise when the observed count is 0, so may be attributed to the numerical instability which arises in the tails.
#include "memory.h" #include <math.h> #include "random.h" xy gaussdisc(double dd,double rr,double tol,int maxeval) ; static double euclid2(xy A,xy B) { return (A.x-B.x)*(A.x-B.x) + (A.y-B.y)*(A.y-B.y) ; } // the proby of getting <n successes from N trials static double ltlbin(int n,int N,double p) { int i ; double u,v,q=1-p,r,s ; for(p/=q,s=r=1,i=1;i<n;i++) s += ( r *= p*(N+1-i)/i ) ; return s*pow(q,N) ; } int main() { int testno,ntests=100000,count,i,n=100000,sno,ex ; xy X ; double q,r,d,qcond,qint ; for(testno=0;testno<ntests;testno++) // if(testno+1==274368) { ranset(testno) ; if((testno+1)%10000==0) printf("%d...\n",testno+1) ; r = 1 / ranf() ; d = 1 / ranf() ; X = gaussdisc(d,r,1e-6,20000) ; qint = X.x ; qcond = X.y ; for(count=sno=0;sno<n;sno++) { X = xy(gaussv(),gaussv()) ; if(euclid2(X,xy(0,d))<r*r) count += 1 ; } if(qint==0||qint==1) { ex = qint?count:0 ; if(count!=ex||qcond>1) printf("%d (cond=%.2f): obs/exp = %d/%d\n",testno,X.y,count,ex) ; } else { q = (count-n*qint)/sqrt(n*qint*(1-qint)) ; if(fabs(q)>4||qcond>1) { if(qint>0&&n*qint<100&&count>n*qint) { q = 1 - ltlbin(count,n,qint) ; if(q<0.00003||qcond>1) printf("%d (cond=%.2f): proby = %.1e%s, obs/exp=%d/%.3f\n", testno,qcond,q,q<3e-7?" ***":"",count,n*qint) ; } else printf("%d (cond=%.2f): sigmage =%4.1f%s, obs/exp=%d/%d\n", testno,qcond,q,fabs(q)>5?" ***":"",count,(int)(0.5+n*qint)) ; } } } }
The results show the similar properties to those of gaussinttest.
gaussint.c gaussinttest.c gaussdisctest.c utils: quadrate.c rtlnorm.c log1plus.c quadrate.h memory.h ranf.c random.h hillclimbing: linmax.c
The following call lines compile and execute the tests, using g++ as compiler:
g++ -o gaussinttest -O -g gaussinttest.c gaussint.c quadrate.c ranf.c rtlnorm.c log1plus.c linmax.c ; gaussinttest g++ -o gaussdisctest -O -g gaussdisctest.c gaussint.c quadrate.c ranf.c rtlnorm.c log1plus.c linmax.c ; gaussdisctest