summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/asymptote/stats.asy
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2013-04-07 18:19:31 +0000
committerKarl Berry <karl@freefriends.org>2013-04-07 18:19:31 +0000
commit752012c605d34cd943795527a9738475a6958fcc (patch)
tree4ee06acdd8333a662c2d6f6ef716235053468f55 /Master/texmf-dist/asymptote/stats.asy
parent9789d09132f18a838e84f041b4b3aff28d3426ec (diff)
texmf -> texmf-dist: start with unique dirs from texmf
git-svn-id: svn://tug.org/texlive/trunk@29712 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/asymptote/stats.asy')
-rw-r--r--Master/texmf-dist/asymptote/stats.asy292
1 files changed, 292 insertions, 0 deletions
diff --git a/Master/texmf-dist/asymptote/stats.asy b/Master/texmf-dist/asymptote/stats.asy
new file mode 100644
index 00000000000..be9efdfa429
--- /dev/null
+++ b/Master/texmf-dist/asymptote/stats.asy
@@ -0,0 +1,292 @@
+private import graph;
+
+real legendmarkersize=2mm;
+
+real mean(real A[])
+{
+ return sum(A)/A.length;
+}
+
+// unbiased estimate
+real variance(real A[])
+{
+ return sum((A-mean(A))^2)/(A.length-1);
+}
+
+real variancebiased(real A[])
+{
+ return sum((A-mean(A))^2)/A.length;
+}
+
+// unbiased estimate
+real stdev(real A[])
+{
+ return sqrt(variance(A));
+}
+
+real rms(real A[])
+{
+ return sqrt(sum(A^2)/A.length);
+}
+
+real skewness(real A[])
+{
+ real[] diff=A-mean(A);
+ return sum(diff^3)/sqrt(sum(diff^2)^3/A.length);
+}
+
+real kurtosis(real A[])
+{
+ real[] diff=A-mean(A);
+ return sum(diff^4)/sum(diff^2)^2*A.length;
+}
+
+real kurtosisexcess(real A[])
+{
+ return kurtosis(A)-3;
+}
+
+real Gaussian(real x, real sigma)
+{
+ static real sqrt2pi=sqrt(2pi);
+ return exp(-0.5*(x/sigma)^2)/(sigma*sqrt2pi);
+}
+
+real Gaussian(real x)
+{
+ static real invsqrt2pi=1/sqrt(2pi);
+ return exp(-0.5*x^2)*invsqrt2pi;
+}
+
+// Return frequency count of data in [bins[i],bins[i+1]) for i=0,...,n-1.
+int[] frequency(real[] data, real[] bins)
+{
+ int n=bins.length-1;
+ int[] freq=new int[n];
+ for(int i=0; i < n; ++i)
+ freq[i]=sum(bins[i] <= data & data < bins[i+1]);
+ return freq;
+}
+
+// Return frequency count in n uniform bins from a to b
+// (faster than the above more general algorithm).
+int[] frequency(real[] data, real a, real b, int n)
+{
+ int[] freq=sequence(new int(int x) {return 0;},n);
+ real h=n/(b-a);
+ for(int i=0; i < data.length; ++i) {
+ int I=Floor((data[i]-a)*h);
+ if(I >= 0 && I < n)
+ ++freq[I];
+ }
+ return freq;
+}
+
+// Return frequency count in [xbins[i],xbins[i+1]) and [ybins[j],ybins[j+1]).
+int[][] frequency(real[] x, real[] y, real[] xbins, real[] ybins)
+{
+ int n=xbins.length-1;
+ int m=ybins.length-1;
+ int[][] freq=new int[n][m];
+ bool[][] inybin=new bool[m][y.length];
+ for(int j=0; j < m; ++j)
+ inybin[j]=ybins[j] <= y & y < ybins[j+1];
+ for(int i=0; i < n; ++i) {
+ bool[] inxbini=xbins[i] <= x & x < xbins[i+1];
+ int[] freqi=freq[i];
+ for(int j=0; j < m; ++j)
+ freqi[j]=sum(inxbini & inybin[j]);
+ }
+ return freq;
+}
+
+// Return frequency count in nx by ny uniform bins in box(a,b).
+int[][] frequency(real[] x, real[] y, pair a, pair b, int nx, int ny=nx)
+{
+ int[][] freq=new int[nx][];
+ for(int i=0; i < nx; ++i)
+ freq[i]=sequence(new int(int x) {return 0;},ny);
+ real hx=nx/(b.x-a.x);
+ real hy=ny/(b.y-a.y);
+ real ax=a.x;
+ real ay=a.y;
+ for(int i=0; i < x.length; ++i) {
+ int I=Floor((x[i]-ax)*hx);
+ int J=Floor((y[i]-ay)*hy);
+ if(I >= 0 && I <= nx && J >= 0 && J <= ny)
+ ++freq[I][J];
+ }
+ return freq;
+}
+
+int[][] frequency(pair[] z, pair a, pair b, int nx, int ny=nx)
+{
+ int[][] freq=new int[nx][];
+ for(int i=0; i < nx; ++i)
+ freq[i]=sequence(new int(int x) {return 0;},ny);
+ real hx=nx/(b.x-a.x);
+ real hy=ny/(b.y-a.y);
+ real ax=a.x;
+ real ay=a.y;
+ for(int i=0; i < z.length; ++i) {
+ int I=Floor((z[i].x-ax)*hx);
+ int J=Floor((z[i].y-ay)*hy);
+ if(I >= 0 && I < nx && J >= 0 && J < ny)
+ ++freq[I][J];
+ }
+ return freq;
+}
+
+path halfbox(pair a, pair b)
+{
+ return a--(a.x,b.y)--b;
+}
+
+path topbox(pair a, pair b)
+{
+ return a--(a.x,b.y)--b--(b.x,a.y);
+}
+
+// Draw a histogram for bin boundaries bin[n+1] of frequency data in count[n].
+void histogram(picture pic=currentpicture, real[] bins, real[] count,
+ real low=-infinity,
+ pen fillpen=nullpen, pen drawpen=nullpen, bool bars=false,
+ Label legend="", real markersize=legendmarkersize)
+{
+ if((fillpen == nullpen || bars == true) && drawpen == nullpen)
+ drawpen=currentpen;
+ bool[] valid=count > 0;
+ real m=min(valid ? count : null);
+ real M=max(valid ? count : null);
+ bounds my=autoscale(pic.scale.y.scale.T(m),pic.scale.y.T(M),
+ pic.scale.y.scale);
+ if(low == -infinity) low=pic.scale.y.scale.Tinv(my.min);
+ real last=low;
+ int n=count.length;
+ begingroup(pic);
+ for(int i=0; i < n; ++i) {
+ if(valid[i]) {
+ real c=count[i];
+ pair b=Scale(pic,(bins[i+1],c));
+ pair a=Scale(pic,(bins[i],low));
+ if(fillpen != nullpen) {
+ fill(pic,box(a,b),fillpen);
+ if(!bars) draw(pic,b--(b.x,a.y),fillpen);
+ }
+ if(!bars)
+ draw(pic,halfbox(Scale(pic,(bins[i],last)),b),drawpen);
+ else draw(pic,topbox(a,b),drawpen);
+ last=c;
+ } else {
+ if(!bars && last != low) {
+ draw(pic,Scale(pic,(bins[i],last))--Scale(pic,(bins[i],low)),drawpen);
+ last=low;
+ }
+ }
+ }
+ if(!bars && last != low)
+ draw(pic,Scale(pic,(bins[n],last))--Scale(pic,(bins[n],low)),drawpen);
+ endgroup(pic);
+
+ if(legend.s != "") {
+ marker m=marker(scale(markersize)*shift((-0.5,-0.5))*unitsquare,
+ drawpen,fillpen == nullpen ? Draw :
+ (drawpen == nullpen ? Fill(fillpen) : FillDraw(fillpen)));
+ legend.p(drawpen);
+ pic.legend.push(Legend(legend.s,legend.p,invisible,m.f));
+ }
+}
+
+// Draw a histogram for data in n uniform bins between a and b
+// (optionally normalized).
+void histogram(picture pic=currentpicture, real[] data, real a, real b, int n,
+ bool normalize=false, real low=-infinity,
+ pen fillpen=nullpen, pen drawpen=nullpen, bool bars=false,
+ Label legend="", real markersize=legendmarkersize)
+{
+ real dx=(b-a)/n;
+ real[] freq=frequency(data,a,b,n);
+ if(normalize) freq /= dx*sum(freq);
+ histogram(pic,a+sequence(n+1)*dx,freq,low,fillpen,drawpen,bars,legend,
+ markersize);
+}
+
+// Method of Shimazaki and Shinomoto for selecting the optimal number of bins.
+// Shimazaki H. and Shinomoto S., A method for selecting the bin size of a
+// time histogram, Neural Computation (2007), Vol. 19(6), 1503-1527.
+// cf. http://www.ton.scphys.kyoto-u.ac.jp/~hideaki/res/histogram.html
+int bins(real[] data, int max=100)
+{
+ real m=min(data);
+ real M=max(data)*(1+epsilon);
+ real n=data.length;
+ int bins=1;
+ real minC=2n-n^2; // Cost function for N=1.
+ for(int N=2; N <= max; ++N) {
+ real C=N*(2n-sum(frequency(data,m,M,N)^2));
+ if(C < minC) {
+ minC=C;
+ bins=N;
+ }
+ }
+
+ return bins;
+}
+
+// return a pair of central Gaussian random numbers with unit variance
+pair Gaussrandpair()
+{
+ real r2,v1,v2;
+ do {
+ v1=2.0*unitrand()-1.0;
+ v2=2.0*unitrand()-1.0;
+ r2=v1*v1+v2*v2;
+ } while(r2 >= 1.0 || r2 == 0.0);
+ return (v1,v2)*sqrt(-log(r2)/r2);
+}
+
+// return a central Gaussian random number with unit variance
+real Gaussrand()
+{
+ static real sqrt2=sqrt(2.0);
+ static pair z;
+ static bool cached=true;
+ cached=!cached;
+ if(cached) return sqrt2*z.y;
+ z=Gaussrandpair();
+ return sqrt2*z.x;
+}
+
+struct linefit {
+ real m,b; // slope, intercept
+ real dm,db; // standard error in slope, intercept
+ real r; // correlation coefficient
+ real fit(real x) {
+ return m*x+b;
+ }
+}
+
+// Do a least-squares fit of data in real arrays x and y to the line y=m*x+b
+linefit leastsquares(real[] x, real[] y)
+{
+ linefit L;
+ int n=x.length;
+ if(n == 1) abort("Least squares fit requires at least 2 data points");
+ real sx=sum(x);
+ real sy=sum(y);
+ real sxx=n*sum(x^2)-sx^2;
+ real sxy=n*sum(x*y)-sx*sy;
+ L.m=sxy/sxx;
+ L.b=(sy-L.m*sx)/n;
+ if(n > 2) {
+ real syy=n*sum(y^2)-sy^2;
+ if(sxx == 0 || syy == 0) return L;
+ L.r=sxy/sqrt(sxx*syy);
+ real arg=syy-sxy^2/sxx;
+ if(arg <= 0) return L;
+ real s=sqrt(arg/(n-2));
+ L.dm=s*sqrt(1/sxx);
+ L.db=s*sqrt(1+sx^2/sxx)/n;
+ }
+ return L;
+}