: \$Id: stats.mod,v 1.227 2012/09/07 14:25:41 samn Exp \$ :* COMMENT COMMENT randwd randomly chooses n bits to set to 1 hamming v.hamming(v1) is hamming distance between 2 vecs flipbits v.flipbits(scratch,num) flips num rand chosen bits flipbalbits v.flipbalbits(scratch,num) balanced flipping vpr v.vpr prints out vector as 1 (x[i]>0) or 0 (x[i]<=0) fac not vec related - returns factorial logfac not vec related - returns log factorial vseed set some C level randomizer seeds slope(num) does a linear regression to find the slope, assuming num=timestep of vector vslope(v2) does a linear regression to find the slope, assuming num=timestep of vector stats(num,[out]) does a linear regression, assuming num=timestep of vector vstats(v2,[out]) does a linear regression, using v2 as the x-coords setrnd(v,flag) does set rand using 1:rand, 2:drand48 v.hash(list) // make a hash out values in vecs in list v.unnan([nan_value,][inf_value]) // remove nan's and inf's from a vector ENDCOMMENT NEURON { SUFFIX stats GLOBAL INSTALLED,seed,kmeasure,verbose,self_ok_combi,hretval,flag,transpose,newline } PARAMETER { : BVBASE = 0. : defined in vecst.mod INSTALLED=0 kmeasure=0 verbose=0 self_ok_combi=0 hretval=0 transpose=0 newline=90 flag=0 : flag can be used by any of the routines for different things } ASSIGNED { seed } VERBATIM #include "misc.h" #define MIN_MERGESORT_LIST_SIZE 32 union dblint { int i[2]; double d; }; static u_int32_t ilow=0; static u_int32_t ihigh=0; // same as ihigh in (/usr/site/nrniv/nrn/src/ivoc/ivocrand.cpp:75) static double *x1x, *y1y, *z1z; static void vprpr(); static int compare_ul(const void* l1, const void* l2) { int retval; unsigned long d; d = (*((const unsigned long*) l1)) - (*((const unsigned long*) l2)); if(d==0) return 1; if(d < 0) return -1; return 0; } ENDVERBATIM :* v1.slope(num) does a linear regression to find the slope, assuming num=timestep of vector VERBATIM static double slope(void* vv) { int i, n; double *x, *y; double timestep, sigxy, sigx, sigy, sigx2; /* how to get the instance data */ n = vector_instance_px(vv, &y); if(ifarg(1)) { timestep = *getarg(1); } else { printf("You must supply a timestep\n"); return 0; } sigxy= sigx= sigy= sigx2=0; // initialize these x = (double *) malloc(sizeof(double)*n); for(i=0; i*maxsqerr) *maxsqerr = val; *meansqerr += val; } *meansqerr=*meansqerr/(double)n; return *meansqerr; } ENDVERBATIM :* v1.stats(num) does a linear regression, assuming num=timestep of vector VERBATIM static double stats(void* vv) { int i, n; double *x, *y, *out; double timestep, sigxy, sigx, sigy, sigx2, sigy2; double r, m, b, dmeansqerr,dmaxsqerr; /* how to get the instance data */ n = vector_instance_px(vv, &y); if(ifarg(1)) { timestep = *getarg(1); } else { printf("You must supply a timestep\n"); return 0; } sigxy= sigx= sigy= sigx2=sigy2= 0; // initialize these x = (double *) malloc(sizeof(double)*n); for(i=0; iX[p->pidse[0]]; Y=&p->Y[p->pidse[1]]; X = p->X; Y = p->Y; p->sigx=p->sigy=p->sigxy=p->sigx2=p->sigy2=0.0; for(i=p->pidse[0]; ipidse[1]; i++) { // p->sigxy += *X * *Y; // p->sigx += *X; // p->sigy += *Y; // p->sigx2 += *X * *X; // p->sigy2 += *Y * *Y; // X++; Y++; p->sigxy += X[i] * Y[i]; p->sigx += X[i]; p->sigy += Y[i]; p->sigx2 += X[i] * X[i]; p->sigy2 += Y[i] * Y[i]; } return NULL; } /* v1.pcorrels2(v2) does a Pearson correlation*/ #if defined(t) static double pcorrelsmt(double *x, double* y, int n,int nth) { int i,nperth,idx,rc; double sigxy, sigx, sigy, sigx2, sigy2, ret; pcorst** pp; pthread_t* pth; pthread_attr_t attr; ret=sigxy=sigx=sigy=sigx2=sigy2=0.0; // initialize these nperth = n / nth; //allocate thread args pp = (pcorst**)malloc(sizeof(pcorst*)*nth); idx=0; for(i=0;iX = x; pp[i]->Y = y; pp[i]->pidse[0] = idx; pp[i]->pidse[1] = idx + nperth; idx += nperth; } i--; if(pp[i]->pidse[1] < n || pp[i]->pidse[1] > n) pp[i]->pidse[1] = n; //make sure all values used //allocate thread IDs pth=(pthread_t*)malloc(sizeof(pthread_t)*nth); pthread_attr_init(&attr); pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); //start threads for(i=0;isigx; sigy += pp[i]->sigy; sigxy += pp[i]->sigxy; sigx2 += pp[i]->sigx2; sigy2 += pp[i]->sigy2; } sigxy -= (sigx * sigy) / n; sigx2 -= (sigx * sigx) / n; sigy2 -= (sigy * sigy) / n; if(sigx2 <= 0) goto PCMTFREE; if(sigy2 <= 0) goto PCMTFREE; ret = sigxy / sqrt(sigx2*sigy2); PCMTFREE: //free memory for(i=0;i 1.0 ){ printf("ppval ERRA: r=%g must be : -1.0 <= r <= 1.0\n",r); return -1.0; } if( n < 3 ){ printf("ppval ERRB: n too small, can't calc probability on samples with < 3 values!\n"); return -1.0; } df = n-2; // degres of freedom // Use a small floating point value to prevent divide-by-zero nonsense // fixme: TINY is probably not the right value and this is probably not // the way to be robust. The scheme used in spearmanr is probably better. TINY = 1.0e-20; ts = r*sqrt(df/((1.0-r+TINY)*(1.0+r+TINY))); mpval = betai(0.5*df,0.5,df/(df+ts*ts)); return mpval; ENDVERBATIM } VERBATIM static const double* sortdata = NULL; /* used in the quicksort algorithm */ /* Helper function for sort. Previously, this was a nested function under * sort, which is not allowed under ANSI C. */ static int compare(const void* a, const void* b) { const int i1 = *(const int*)a; const int i2 = *(const int*)b; const double term1 = sortdata[i1]; const double term2 = sortdata[i2]; if (term1 < term2) return -1; if (term1 > term2) return +1; return 0; } void csort (int n, const double mdata[], int index[]) /* Sets up an index table given the data, such that mdata[index[]] is in * increasing order. Sorting is done on the indices; the array mdata * is unchanged. */ { int i; sortdata = mdata; for (i = 0; i < n; i++) index[i] = i; qsort(index, n, sizeof(int), compare); } static double* getrank (int n, double mdata[]) /* Calculates the ranks of the elements in the array mdata. Two elements with * the same value get the same rank, equal to the average of the ranks had the * elements different values. The ranks are returned as a newly allocated * array that should be freed by the calling routine. If getrank fails due to * a memory allocation error, it returns NULL. */ { int i; double* rank; int* index; rank = malloc(n*sizeof(double)); if (!rank) return NULL; index = malloc(n*sizeof(int)); if (!index) { free(rank); return NULL; } /* Call csort to get an index table */ csort (n, mdata, index); /* Build a rank table */ for (i = 0; i < n; i++) rank[index[i]] = i; /* Fix for equal ranks */ i = 0; while (i < n) { int m; double value = mdata[index[i]]; int j = i + 1; while (j < n && mdata[index[j]] == value) j++; m = j - i; /* number of equal ranks found */ value = rank[index[i]] + (m-1)/2.; for (j = i; j < i + m; j++) rank[index[j]] = value; i += m; } free (index); return rank; } /* The spearman routine calculates the Spearman rank correlation between two vectors. n (input) int The number of elements in a data vector data1 (input) double array -- the first vector data2 (input) double array -- the second vector */ static double spearman(int n, double* data1, double* data2) { int i; int m = 0; double* rank1; double* rank2; double result = 0.; double denom1 = 0.; double denom2 = 0.; double avgrank; double* tdata1; double* tdata2; tdata1 = malloc(n*sizeof(double)); if(!tdata1) return 0.0; /* Memory allocation error */ tdata2 = malloc(n*sizeof(double)); if(!tdata2) /* Memory allocation error */ { free(tdata1); return 0.0; } for (i = 0; i < n; i++) { tdata1[m] = data1[i]; tdata2[m] = data2[i]; m++; } if (m==0) return 0; rank1 = getrank(m, tdata1); free(tdata1); if(!rank1) return 0.0; /* Memory allocation error */ rank2 = getrank(m, tdata2); free(tdata2); if(!rank2) /* Memory allocation error */ { free(rank1); return 0.0; } avgrank = 0.5*(m-1); /* Average rank */ for (i = 0; i < m; i++) { const double value1 = rank1[i]; const double value2 = rank2[i]; result += value1 * value2; denom1 += value1 * value1; denom2 += value2 * value2; } /* Note: denom1 and denom2 cannot be calculated directly from the number * of elements. If two elements have the same rank, the squared sum of * their ranks will change. */ free(rank1); free(rank2); result /= m; denom1 /= m; denom2 /= m; result -= avgrank * avgrank; denom1 -= avgrank * avgrank; denom2 -= avgrank * avgrank; if (denom1 <= 0) return 0; /* include '<' to deal with roundoff errors */ if (denom2 <= 0) return 0; /* include '<' to deal with roundoff errors */ result = result / sqrt(denom1*denom2); return result; } static double scorrel(void* vv) { int i, n; double *x, *y; n = vector_instance_px(vv, &x); if ((i=vector_arg_px(1, &y)) != n ) {printf("scorrelERRA: %d %d\n",n,i); hxe();} return spearman(n,x,y); } //* Kendall's correlation routines //** Erfcc() Returns the complementary error function erfc(x) with fractional error // everywhere less than 1.2 x 10^-7. // from place.mod, is that from numerical recipes?? double Erfcc (double x) { double mt,z,ans; z=fabs(x); mt=1.0/(1.0+0.5*z); ans=mt*exp(-z*z-1.26551223+mt*(1.00002368+mt*(0.37409196+mt*(0.09678418+\ mt*(-0.18628806+mt*(0.27886807+mt*(-1.13520398+mt*(1.48851587+\ mt*(-0.82215223+mt*0.17087277))))))))); return x >= 0.0 ? ans : 2.0-ans; } //** Rktau() R version of kendall tau, doesnt have huge memory footprint //function returns kendall's tau double Rktau (double* x, double* y, int n){ int i,j; double c,vx,vy,sx,sy,var,z,tau; c = vx = vy = 0.0; for(i = 0; i < n; i++) { for(j = 0; j < i; j++) { sx = (x[i] - x[j]); sx = ((sx > 0) ? 1 : ((sx == 0)? 0 : -1)); sy = (y[i] - y[j]); sy = ((sy > 0) ? 1 : ((sy == 0)? 0 : -1)); vx += sx * sx; vy += sy * sy; c += sx * sy; } } if(vx>0 && vy>0) { tau = c / sqrt(vx*vy); return tau; } return 0.; } //** vec1.kcorrel(vec2,[fast version -- useful for large arrays, vec of size 1 holding p-value]) // kendall's tau correlation static double kcorrel (void* vv) { int i, n; double *x, *y, *prob, *i1d, *i2d, *ps, var, z, tau; n = vector_instance_px(vv, &x); if ((i=vector_arg_px(1, &y)) != n ) {printf("kcorrel ERRA: %d %d\n",n,i); hxe();} if(ifarg(2) && *getarg(2)) { i1d=dcrset(n*3); i2d=&i1d[n]; ps=&i2d[n]; tau=kcorfast(x,y,i1d,i2d,n,ps); } else { tau = Rktau(x,y,n); } if(!(ifarg(3) && vector_arg_px(3,&prob))) prob = 0x0; //does user want to store p-value? if(prob) { //get p-value var = (4.0 * n + 10.0) / (9.0 * n * (n - 1.0)); z = tau / sqrt(var); *prob = Erfcc(fabs(z)/1.4142136); //when prob small, chance of tau having its value by chance is small } return tau; } //** mycompare() comparison function for qsort -- sorts in ascending order static int mycompare (const void* a, const void* b) { double d1,d2; d1 = *(double*)a; d2 = *(double*)b; if (d1 < d2) return -1; if (d1 > d2) return +1; return 0; } //** mergesort_array() // recursive mergesort -- sorts a by splitting into sublists, sorting, and recombining void mergesort_array (double a[], int size, double temp[],unsigned long* swapcount) { int i1, i2, i, tempi, j, vv; double *right,*left; if(size<=1) return; //base case -- 1 element is sorted by definition if (0 && size = 0; j--) { if (a[j] <= vv) break; a[j + 1] = a[j]; } a[j + 1] = vv; } return; } mergesort_array(a, size/2, temp,swapcount); //sort left half mergesort_array(a + size/2, size - size/2, temp,swapcount); //sort right half //merge halves together i=tempi=0; i1 = size/2; i2 = size - size/2; left = a; right = &a[size/2]; while(i1>0 && i2>0) { if(*right < *left) { *swapcount += i1; temp[i] = *right++; i2--; } else { temp[i] = *left++; i1--; } i++; } if(i2>0) { while(i2-->=0 && i=0 && i2) printf("nPair=%lu\n",nPair); if(verbose>3){printf("i1d after qsort2: "); for(i=0;i 0) { qsort(&i2d[i-tieCount-1],tieCount+1,sizeof(double),mycompare); m1 += tieCount * (tieCount + 1) / 2; s += getMs(&i2d[i-tieCount-1],tieCount+1); tieCount = 0; } } if(verbose>2) printf("tieCount=%lu\n",tieCount); if(tieCount > 0) { qsort(&i2d[n-tieCount-1],tieCount+1,sizeof(double),mycompare); m1 += tieCount * (tieCount + 1) / 2; s += getMs(&i2d[n-tieCount-1],tieCount+1); } if(verbose>2) printf("tieCount=%lu\n",tieCount); swapCount = 0; mergesort_array(i2d,n,ps,&swapCount); //sort input2 & count # of swaps to get into sorted order if(verbose>3) { printf("i2d after mergesort: "); for(i=0;i2) printf("swapCount=%lu\n",swapCount); m2 = getMs(i2d,n); if(verbose>2) printf("s=%lu m1=%lu m2=%lu\n",s,m1,m2); s -= (m1 + m2) + 2 * swapCount; denom1=nPair-m1; denom2=nPair-m2; if(verbose>2) printf("s=%lu d1=%g d2=%g\n",s,denom1,denom2); if(denom1>0. && denom2>0.) return s / sqrt(denom1*denom2); else return 0.; } //root mean square of vector's elements static double rms (void* vv) { int i,n; double *x,sum; if(!(n=vector_instance_px(vv, &x))) {printf("rms ERRA: 0 sized vector!\n"); hxe();} sum=0.0; for(i=0;i0.) return sqrt(sum); else return 0.0; } //cumulative sum of vector's elements static double cumsum (void* vv) { int i,n; double *x,*y; if(!(n=vector_instance_px(vv, &x))) {printf("cumsum ERRA: 0 sized vector!\n"); hxe();} if(vector_arg_px(1, &y) != n) {printf("cumsum ERRB: output vec size needs size of %d\n",n); hxe();} memcpy(y,x,sizeof(double)*n); for(i=1;i=nx) { hoc_execerror("# of flips exceeds (or ==) vector size", 0); } for (i=0; i < nx; i++) { x[i] = BVBASE; } for (i=0,jj=0; i < flip; i++) { /* flip these bits */ ii = (int) ((nx+1)*drand48()); if (x[ii]==BVBASE) { x[ii] = 1.; if (flag) { y[jj] = ii; jj++; } } else { i--; } } return flip; } ENDVERBATIM :* v1.hash(veclist) VERBATIM static double hash (void* vv) { int i, j, nx, nv[VRRY], num, vfl; union dblint xx; Object* ob; double *x, *vvo[VRRY], big, y, prod; nx = vector_instance_px(vv, &x); if (ifarg(1)) { vfl=0; ob=*hoc_objgetarg(1); num = ivoc_list_count(ob); if (num>VRRY) {printf("vecst:hash ERR: can only handle %d vecs: %d\n",VRRY,num); hxe();} for (i=0;iVRRY) {printf("vecst:smash ERRA: can only handle %d vecs: %d\n",VRRY,num); hxe();} if (transpose) if (nx!=num) { printf("vecst:smash ERRB %d %d %d\n",i,nx,nv[i]);hxe(); } for (i=0;iVRRY) {printf("stats:dpro ERR: can only handle %d vecs: %d\n",VRRY,num); hxe();} for (i=0;i= randvar and mul by step : v1.setrnd(4.5,min,max[,seed]) // [min,max) : v1.setrnd(5[,n,seed]) -- integers [0,100) or [0,n] : v1.setrnd(5[,min,max,seed]) -- integers [min,max] -- if seed=0 it's not reset : v1.setrnd(5,ind[,seed]) -- random values from ind : v1.setrnd(6) -- unique integers as follows: : v1.setrnd(6,min,max[,seed]) -- unique in [min,max] : v1.setrnd(6,min,max,exclude_vec[,seed]) -- unique in [min,max] excluding values in exclude_vec VERBATIM static double setrnd (void* vv) { int flag, i,j,k,n,cnt; unsigned int nx, nx1, nex, lt, rt, mid; double *x, y, *ex, *ex2, min, max, dfl, tmp, step, num; unsigned long value; value=1; nx = vector_instance_px(vv, &x); flag = (int)(dfl=*getarg(1)); if (flag==1) { for (i=0; i < nx; i++) x[i] = (double)rand()/RAND_MAX; } else if (flag==2) { for (i=0; i < nx; i++) x[i] = drand48(); } else if (flag==3) { // scop_random()'s cheap and dirty rand unsigned long a = 2147437301, c = 453816981, m = ~0; value = (unsigned long) seed; for (i=0; i < nx; i++) { value = a * value + c; x[i] = (fabs((double) value / (double) m)); } seed=(double)value; } else if (flag==4) { // mcell_ran4() doubles ex=0x0; i=2; if (ifarg(i)) { if (hoc_is_object_arg(i)) { nex=vector_arg_px(i++,&ex); // vector to look in step=ifarg(i)?*getarg(i):1.0; max=1.0; i++; } else { if (dfl==4.5 || ifarg(4)) { // flag 4.5 resolves ambiguity of arg3 max or seed min=*getarg(i++); max=*getarg(i++)-min; dfl=4.5; } else { max=*getarg(i++); } } } else max=1.0; // default if (ifarg(i)) { y=*getarg(i++); if (y) ihigh=(unsigned int)y; } // look for seed if (max==0) { for (i=0;iex[mid]) { if (num1) { // sort the exclude vector scrset(nex); x1x = (double *)realloc(x1x,sizeof(double)*nx*4); for (i=0;imax || ex[j]0 && ex[j]<=ex[j-1])) { printf("%g in exclusion list -- out of range [%g,%g] or a repeat\n",ex[j],min,max);hxe();} if (max-min+1-nex==nx) { for (i=min,k=0;i<=max;i++) { y=(double)i; for (j=0;j= to an excluded value must shift by # its >= to for (i=0;i=ex[j]) k++; // will move it up by k x[i]=z1z[i]+k; } } else for (i=0;i1;) { mcell_ran4(&ihigh, y, 1, n); n--; k=(int)y[0]; // random int(n) // 0 <= k < n. temp = x[n]; x[n] = x[k]; x[k] = temp; } } //shuffle array of unsigned ints void uishuffle(unsigned int* x,int nx) { int n,k; unsigned int temp; double y[1]; for (n=nx;n>1;) { mcell_ran4(&ihigh, y, 1, n); n--; k=(int)y[0]; // random int(n) // 0 <= k < n. temp = x[n]; x[n] = x[k]; x[k] = temp; } } //shuffle array of unsigned ints void ishuffle(int* x,int nx) { int n,k,temp; double y[1]; for (n=nx;n>1;) { mcell_ran4(&ihigh, y, 1, n); n--; k=(int)y[0]; // random int(n) // 0 <= k < n. temp = x[n]; x[n] = x[k]; x[k] = temp; } } unsigned long choose (int n, int k) { int i,delta; unsigned long ret; //assert ((n >= 0) && (k >= 0)); if (n < k) return 0; if (n==k) return 1; if (k < n - k) { delta = n - k; } else { delta = k; k = n - k; } ret = delta + 1; for (i = 2; i <= k; ++i) ret = (ret * (delta + i)) / i; return ret; } // computes combinadic given combination vector unsigned long syncci (int nn,int kk,int* ccvv) { unsigned long c = 0; while ((kk > 0) && (*ccvv >= kk)) { c += choose (*ccvv++, kk--); } return c; } // computes combinadic given combination vector unsigned long syncc (int nn,int kk,double* ccvv) { unsigned long c = 0; while ((kk > 0) && (*ccvv >= kk)) { c += choose (*ccvv++, kk--); } return c; } // computes combination vector given combinadic void synccv (int nn,int kk,int cc,double* ccvv) { unsigned long n_k; while (--nn >= 0) { n_k = choose (nn, kk); if (cc >= n_k) { cc -= n_k; *ccvv++ = nn; --kk; } } } ENDVERBATIM : returns the # of combinations for choosing k from a set of n : combs(n,k) FUNCTION combs () { VERBATIM unsigned int n,k; n=(unsigned int)*getarg(1); k=(unsigned int)*getarg(2); return choose(n,k); ENDVERBATIM } :* vec.comb(nn,kk,cc) : returns combination indexed by cc, from selecting kk elements from set of nn elements : element ids are from 0..nn-1 VERBATIM static double comb (void* vv) { double* x; int nn,kk,cc,sz,i; sz=vector_instance_px(vv,&x); nn=(int)*getarg(1); kk=(int)*getarg(2); cc=(int)*getarg(3); if(sz= %d , is %d!\n",kk,sz); return 0.0; } memset(x,0,sizeof(double)*kk); synccv(nn,kk,cc,x); vector_resize(vv,kk); return 1.0; } ENDVERBATIM :* vec.combid(nn,kk) : returns combination index corresponding to data in vector : indices from selecting kk elements from set of nn elements : element ids are from 0..nn-1 VERBATIM static double combid (void* vv) { double* x; int nn,kk,sz,i; sz=vector_instance_px(vv,&x); nn=(int)*getarg(1); kk=(int)*getarg(2); if(sz= %d , is %d!\n",kk,sz); return 0.0; } return syncc(nn,kk,x); } ENDVERBATIM VERBATIM int findlong (unsigned long* p,unsigned long val,int istart,int iend) { int i; for(i=istart;i<=iend;i++) if(p[i]==val) return 1; return 0; } ENDVERBATIM :* rsampsig(vIN0,vIN1,prc,"vhfmeasure","compfunc",vhsout,[onesided,nocmbchk,inorder]) : vIN0 - elements of group 0 : vIN1 - elements of group 1 : prc - fraction of subsets to measure with "vhfmeasure" , or total # of combinations to test iff prc > 1.0 : "vhfmeasure" - hoc function that takes 1 vector as arg and returns a double, must set hretval_stats : to the return value so can be read from here : "compfunc" - hoc function that takes 2 doubles as args & returns 1 iff a boolean condition is satisfied, must set : hretval_stats to return value so can be read from here : vhsout - vector used in hocmeasure , must have size >= vIN0.size + vIN1.size : onesided - optional, use one sided p-value or two-sided (default = 1) , onsided==0 means use two-sided : nocmbchk - skip combination ID check, useful for groups with large size ( > 30000 ) , to avoid overflow, def=1 == skip combination ID check : returns -1.0 on error VERBATIM static double rsampsig(void* vv){ int n0,n1,na,nn,kk,cc,i,j,*pm,szthis,onesided,nocmbchk,bti,*pids; unsigned long nruncombs,nallcombs,*pcombids; double *g0,*g1,*ga,prc,*g0t,*g1t,dmobs,dm0,dm1,*phso,nmatch,*pthis,dret; void* vhso; //vector * for changing size at end Symbol* pHocVecFunc,*pHocCompFunc; //hoc function pointers dret=-1.0; g0t=g1t=NULL; pm=pids=NULL; pcombids=NULL;//init arrays to null szthis=vector_instance_px(vv, &pthis); //size of calling vector n0=vector_arg_px(1,&g0);//group 0 size n1=vector_arg_px(2,&g1);//group 1 size na=n0+n1;//total # of elements prc=*getarg(3);//% of combinations to try, or total # of combinations to try iff > 1.0 if(!(pHocVecFunc=hoc_lookup(gargstr(4)))){ printf("rsampsig ERRA: couldn't find hoc vec func %s\n",gargstr(4)); goto CLEANRSAMPSIG; } if(!(pHocCompFunc=hoc_lookup(gargstr(5)))){ printf("rsampsig ERRB: couldn't find hoc comp func %s\n",gargstr(5)); goto CLEANRSAMPSIG; } if(vector_arg_px(6,&phso)= %d!\n",na); goto CLEANRSAMPSIG; } if(prc<=0.0) { printf("rsampsig ERRD: invalid value for arg 3, must be > 0.0!\n"); goto CLEANRSAMPSIG; } vhso=vector_arg(6);//get vector arg for resize ga=(double*)malloc(sizeof(double)*na);//all elements memcpy(ga,g0,sizeof(double)*n0);//copy elements of group 0 to ga memcpy(ga+n0,g1,sizeof(double)*n1);//append elements of group 1 to ga //form nruncombs combinations, select elements into groups, and compare measure on groups g0t=(double*)malloc(sizeof(double)*n0);//temp storage for group 0 g1t=(double*)malloc(sizeof(double)*n1);//temp storage for group 1 if(n01) printf("choose(%d,%d)=%ld\n",na,kk,choose(na,kk)); nallcombs=choose(na,kk);//all possible combinations selecting kk elements from a set of size na nruncombs=prc>1.0?prc:prc*nallcombs; nmatch=0.0;//# of combinations to run if(szthis1) printf("na=%d , kk=%d, n0=%d, n1=%d\n",na,kk,n0,n1); if(verbose>1) printf("nruncombs = %ld\n",nruncombs); if(verbose>2) pm=(int*)malloc(sizeof(int)*na); for(i=0;i0 && findlong(pcombids,pcombids[i],0,i-1)); //make sure pcombids[i] wasnt used yet if(verbose) if(i%100==0) { printf("."); fflush(stdout); } for(j=0;j2){ for(j=0;j observed test statistic, to get p-value for(i=0;i dmobs) nmatch++; dret=nmatch/(double)nruncombs;//this output value doesnt mean anything yet CLEANRSAMPSIG: //free memory and return if(ga) free(ga); if(g0t) free(g0t); if(g1t) free(g1t); if(pm) free(pm); if(pcombids) free(pcombids); if(pids) free(pids); return dret; } ENDVERBATIM VERBATIM // rantran(index_vec1,value_vec1[,index_vec2,value_vec2,...]) // index_vec can be replaced with a step which gives regular pointers or base values // in presence of a step value can be eg 10 for 0..9 [or 10.3 for 0,3,6,9 -- not implemented] // rantran() translates a random index vec through a series of mappings from a value onto // one or more possible values: eg rand col indices to rand minicols within that col to // rand cell within the minicol to rand locations on the cell dend static double rantran (void* vv) { int i,j,ix,ixe,ixvn,nvn,rvn,na,xj; double *ixv, *nv, *x, y[1], ixn,step,indx; rvn=vector_instance_px(vv, &x); for (na=1;ifarg(na);na++); na--; // count args for (i=1;i1.) { x1x = (double *)realloc(x1x,sizeof(double)*rvn); mcell_ran4(&ihigh, x1x, rvn, indx); } } for (j=0;j-1) { // step in by xi[type] and then augment by random val x[j]=step + x[j]*indx + ((indx>1.)?(floor(x1x[j])):0.0); } else { xj=(int)x[j]; // prior level index ix=(int)ixv[xj]; ixn=ixv[xj+1]-ix; // possible next level indices; pick 1 in [ix,ixe] if (ixn==1.) { x[j]=nv[ix]; } else { mcell_ran4(&ihigh, y, 1, ixn); // pick 1 rand value [0,ix) x[j]=nv[(int)y[0]+ix]; } } } } return (double)rvn; // number of substitutions performed } ENDVERBATIM :* vec.shuffle() Fisher-Yates shuffle (from wikipedia) : pointlessly extended to augment the vector n-fold (used for cell projections in columns) : makes more sense to just have a parallel random vector of intra-columnar indices VERBATIM static double shuffle (void* vv) { int i,j,k,n,nx,augfac; double *x, y[1], temp, augstep; nx=vector_instance_px(vv, &x); if (ifarg(1)) { augfac=(int)*getarg(1); if (ifarg(2)) augstep=*getarg(2); else augstep=1.0/augfac; x=vector_newsize(vv,nx*augfac); for (i=1;inx) {printf("stats vpr ERRA OOB: %d %d\n",min,max); hxe();} if (ifarg(1)) { if (hoc_is_double_arg(1)) { flag=(int) *getarg(1); if (flag==2) { // binary which is default for (i=min; iBVBASE)?1:0); } else if (flag==0) { for (i=min; i=10?"|":"",(int)x[i],x[i]>=10?"|":""); } else if (flag==-1) { for (i=min; iBVBASE) { fprintf(f,"%d",1); } else { fprintf(f,"%d",0); } } fprintf(f,"\n"); } } else { for (i=min; iBVBASE)?1:0); printf("\n"); } return 1.; } ENDVERBATIM :* v1.vpr2(v2,[BASE]) prints out 2 vectors using = where same -- optional arg allows base 10,16,64 : use VERBOSE to adjust printing verbose=1 for no spaces; verbose=2 no spaces and abbrev 00->_ VERBATIM static double vpr2 (void* vv) { int i, flag, n, nx, ny, colc, base, min,max, fl2; double *x, *y, cnt, ign; char c; nx = vector_instance_px(vv, &x); cnt=min=0; max=nx; ny = vector_arg_px(1, &y); if (nx!=ny) {printf("vpr2 diff sizes %d %d\n",nx,ny); hxe();} base=(int)*getarg(2); flag=(ifarg(3)?(int)*getarg(3):2); for (i=min,fl2=0,colc=0; i0 && x[i]==BVBASE && y[i]==BVBASE) { if (!fl2) {printf(" _ "); colc+=3;} fl2=1; } else { fl2=0; colc++; vprpr(x[i],base); if (colc>(int)newline){printf("\n "); colc=0;} } } printf("\n"); for (i=min,fl2=0,colc=0; i0 && x[i]==BVBASE && y[i]==BVBASE) { if (!fl2) {printf(" _ "); colc+=3;} fl2=1; } else { fl2=0; vprpr(y[i],base); colc++; if (colc>(int)newline){printf("\n ",colc); colc=0;} } } printf("\n"); if (flag==1) { for (i=min,n=0,fl2=0,colc=0; i(int)newline){printf("\n ",colc); colc=0;} } printf("\n"); } } static void vprpr (double x, int base) { int xx; xx=(int)x; if (base==0) { printf("%5.2f",x); } else if (xx>=base && base!=0) { printf("+"); } else if (base==64) { // base 64 if (xx<16) { printf("%x",xx); // 0-f 0-15 } else if (xx<36) { printf("%c",xx+87); // g-z 16-35 } else if (xx<62) { printf("%c",xx+29); // A-Z 36-61 } else if (xx<63) { printf("@"); // @ 62 } else if (xx<64) { printf("="); // = 63 } else printf("ERROR"); } else if (base==10) { printf("%d",xx); } else if (base==16) { printf("%x",xx); } } ENDVERBATIM :* v1.bin(targ,invl[min,max]) place counts for each interval :* v1.bin(list,invl[min,max]) using NQS("count","index") for easy sorting : like .hist() but doesn't throw away values max : note that optional max denotes the start of an interval not end VERBATIM static double bin (void* vv) { int i, j, nx, maxsz, lfl; double* x, *y, *ix, invl, min, max, maxf, jj; Object* ob; void* voi[2]; min=0; max=1e9; maxf=-1e9; nx = vector_instance_px(vv, &x); ob = *hoc_objgetarg(1); if (strncmp(hoc_object_name(ob),"Vector",6)==0) { lfl=0; // lfl is list flag if ((maxsz=openvec(1, &y))==0) hxe(); voi[0]=vector_arg(1); } else { // list of 2 lfl=1; maxsz = list_vector_px3(ob, 0, &y, &voi[0]); if (maxsz!=(i=list_vector_px3(ob,1,&ix,&voi[1]))){printf("binERRA %d %d\n",maxsz,i); hxe();} } invl = *getarg(2); if (ifarg(4)) {min=*getarg(3); max=*getarg(4); } else if (ifarg(3)) max=*getarg(3); for (j=0; j=max) { y[(j=(int)(jj=(max-min)/invl))]++; } else if (x[i]<=min) { y[(j=0)]++; jj=0.; } else { if (x[i]>maxf) maxf=x[i]; // max found j=(int)(jj=(x[i]-min)/invl); if (j>maxsz-1) {printf("bin ERRB OOB: %d>%d\n",j,maxsz-1); hxe();} y[j]++; } if (lfl) ix[j]=jj+min; } maxsz=(max==1e9)?(int)(maxf/invl+1):(int)((max-min)/invl+1); vector_resize(voi[0], maxsz); if (lfl) vector_resize(voi[1], maxsz); return (double)maxsz; } ENDVERBATIM :* ind.ihist(vec,list,tmin,tmax,binsz,[,ioff]) does binning for individual indices : into list -- ioff gives an offset such that index N maps to list.o(0) : normally used for a raster pair of cell indices (ind) and spike time (tvec) : eg nq=new NQS(#cells) : ind.ihist(tvec,nq.vl,tmin,tmax,100) VERBATIM static double ihist (void* vv) { unsigned int i, j, k, n, nx, c; double *x, *tv, ioff, min, max, nbin, binsz; ListVec* pL; Object* obl; nx = vector_instance_px(vv, &x); // vector of indices i = vector_arg_px(1, &tv); // vector of times if (i!=nx) {printf("vecst:ihist()ERR0: diff size %d %d\n",nx,i); hxe();} if (!flag && !ismono1(tv,nx,1)){ printf("vecst:ihist()ERR0A: set flag_stats for non-monotonic time vec\n",nx,i); hxe();} pL = AllocListVec(obl=*hoc_objgetarg(2)); min=*getarg(3); max=*getarg(4); binsz=*getarg(5); if (binsz<=0) {printf("stats:ihist()ERR0B: binsz must be >0 (%g)\n",binsz); hxe();} nbin=floor((max-min)/binsz); max=min+binsz*nbin; // new max if (verbose) printf("%g-%g in %g bins of %g\n",min,max,nbin,binsz); if (ifarg(6)) ioff=*getarg(6); else ioff=0.; if (pL->isz<2) {printf("stats:ihist()ERRA: %d\n",pL->isz); FreeListVec(&pL); hxe();} c=pL->isz; // number of columns (list length) // pL->pv[i][j] -- i goes through the list and j goes through each vector for (i=0;ipv[i]=list_vector_resize(obl, i, (int)nbin); for (j=0;j<(int)nbin;j++) pL->pv[i][j]=0.; } i=n=0; if (!flag) for (;i=min) break; // tvec sorted so can move forward to begin for (;i=max || tv[i]=max) break; k=(int)(x[i]-ioff); // which cell index if (k>=c || k<0) continue; // ignore outside of index range j=(int)((tv[i]-min)/binsz); // which time bin // if(j>=nbin||j<0){printf("INTERR %d %d %d %g %g %g %g\n",k,c,j,nbin,tv[i],min,binsz);hxe();} pL->pv[k][j]++; n++; } flag=0; FreeListVec(&pL); return (double)n; } ENDVERBATIM :* vec.irate(vhist,binsz) - sets contents of vec to instantaneous rate using inverse : of interspike intervals VERBATIM static double irate (void* vv) { unsigned int i, j, n, nx; double *prate,*phist,binsz,t1,t2; nx = vector_arg_px(1, &phist); vector_resize(vv,nx); vector_instance_px(vv, &prate); binsz = *getarg(2); for(i=0;i0) { t1=t2=i; prate[i]=phist[i]*1e3/binsz; i++; break; } } for(;i0) { t2=i; prate[i]=phist[i]*1e3/binsz; break; } } if(verbose>1) if(t1==t2) printf("t1==t2!\n"); for(i=t1;i0) { prate[i] = phist[i]*1e3/binsz; t1=t2; t2=i; if(verbose>2) printf("t1 %g t2 %g\n",t1,t2); } else { if(verbose>1) if(t1==t2) printf("t1==t2!\n"); prate[i] = 1e3 / (binsz*(t2-t1)); } } return 1.0; } //* dlvar(p,sz) -- calculates the local variation in interspike interval stored in p. // p should have all positive numbers. for poisson interspike intervals, dlvar will return 1. // for regular interspike intervals, dlvar will return 0. double dlvar (double* p, int sz) { int i; double s,n,d; s=0.0; for(i=0;i= 2, returning 0!\n"); return 0.0; } return dlvar(x,sz); } ENDVERBATIM VERBATIM //cholesky decomposition from numerical recipies chapter 2.9 //Given a positive-dfinite symmetric matrix a[1..n][1..n], this routine constructs its Cholesky //decomposition, A = L * LT . On input, only the upper triangle of a need be given; it is not //modified. The Cholesky factor L is returned in the lower triangle of a, except for its diagonal //elements which are returned in p[1..n] int choldc(double **a, int n, double p[]) { int i,j,k; double sum; for (i=0;i=0;k--) sum -= a[i][k]*a[j][k]; if (i == j) { if (sum <= 0.0) return 0; p[i]=sqrt(sum); } else a[j][i]=sum/p[i]; } } return 1; } ENDVERBATIM : mktscor(pearson r-value,list of time-series) : each time-series should be normally distributed : upon return, the values will be changed and each pair of time-series will be correlated with pcorrel = ~r : same method as here http://comisef.wikidot.com/tutorial:correlation FUNCTION mktscor () { VERBATIM double dret,*ptmp,**cF,**Y,**X,**cFT,r; ListVec *pTS; int i,j,k,nrows,ncols,nrcf; dret=0.0; pTS=0x0; ptmp=0x0; cF=cFT=0x0; X=Y=0x0; r = *getarg(1); pTS = AllocListVec(*hoc_objgetarg(2)); nrows=pTS->plen[0]; ncols=pTS->isz; X = getdouble2D(nrows,ncols); for(i=0;ipv[i][j]; nrcf=ncols; cF = getdouble2D(nrcf,nrcf); for(i=0;ii) cF[i][j]=0.0; else if(j==i) cF[i][j]=ptmp[i]; cFT = getdouble2D(nrcf,nrcf); for(i=0;ipv[j][i]=Y[i][j]; MKTSCORFREE: dret=1.0; if(pTS) FreeListVec(&pTS); if(ptmp) free(ptmp); if(cF) freedouble2D(&cF,nrcf); if(cFT) freedouble2D(&cFT,nrcf); if(Y) freedouble2D(&Y,nrows); if(X) freedouble2D(&X,nrows); return dret; ENDVERBATIM } :* PROCEDURE install_stats() PROCEDURE install () { if (INSTALLED==1) { printf("\$Id: stats.mod,v 1.227 2012/09/07 14:25:41 samn Exp \$") } else { INSTALLED=1 VERBATIM x1x=y1y=z1z=0x0; ihigh=1; install_vector_method("slope", slope); install_vector_method("moment", moment); install_vector_method("vslope", vslope); install_vector_method("stats", stats); install_vector_method("pcorrel", pcorrel); install_vector_method("scorrel",scorrel); install_vector_method("kcorrel",kcorrel); install_vector_method("rms",rms); install_vector_method("vstats", vstats); install_vector_method("randwd", randwd); install_vector_method("hamming", hamming); install_vector_method("flipbits", flipbits); install_vector_method("flipbalbits", flipbalbits); install_vector_method("vpr", vpr); install_vector_method("vpr2", vpr2); install_vector_method("bin", bin); install_vector_method("ihist", ihist); install_vector_method("setrnd", setrnd); install_vector_method("rantran", rantran); install_vector_method("distance", distance); install_vector_method("ndprd", ndprd); install_vector_method("hash", hash); install_vector_method("smash", smash); install_vector_method("smash1", smash1); install_vector_method("dpro", dpro); install_vector_method("unnan", unnan); install_vector_method("combi", combi); install_vector_method("shuffle", shuffle); install_vector_method("comb", comb); install_vector_method("combid", combid); install_vector_method("rsampsig",rsampsig); install_vector_method("irate",irate); install_vector_method("cumsum",cumsum); install_vector_method("lvar",lvar); ENDVERBATIM } } PROCEDURE prhash (x) { VERBATIM { long unsigned int xx; xx=*(long unsigned int*)(&_lx); printf("%16lx\n",xx); } ENDVERBATIM } :* fac (n) : from numerical recipes p.214 FUNCTION fac (n) { VERBATIM { static int ntop=4; static double a[101]={1.,1.,2.,6.,24.}; static double cof[6]={76.18009173,-86.50532033,24.01409822, -1.231739516,0.120858003e-2,-0.536382e-5}; int j,n; n = (int)_ln; if (n<0) { hoc_execerror("No negative numbers ", 0); } if (n>100) { /* gamma function */ double x,tmp,ser; x = _ln; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.0; for (j=0;j<=5;j++) { x += 1.0; ser += cof[j]/x; } return exp(-tmp+log(2.50662827465*ser)); } else { while (ntop100) { /* gamma function */ double x,tmp,ser; x = _ln; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.0; for (j=0;j<=5;j++) { x += 1.0; ser += cof[j]/x; } return (-tmp+log(2.50662827465*ser)); } else { while (ntopUINT_MAX) y/=1e9; // UINT_MAX is 4.294967e+09 ihigh*=(unsigned int)y; // keep multiplying these out } return ihigh; // global } ENDVERBATIM FUNCTION hashseed () { VERBATIM int i,na,nb,sf=0; double *x=0; if (hoc_is_double_arg(2)) { for (na=1;ifarg(na);na++); na--; if (na==1) nb=2; else nb=na; x = (double *) malloc(sizeof(double)*nb); for (i=1;i<=na;i++) x[i-1]=*getarg(i); if (na==1) x[1]=x[0]*x[0]+13; // so can get something sf=1;//should free } else { nb=(int)*getarg(1); x=hoc_pgetarg(2); } hashseed2(nb,x); if(sf) free(x); //only free if malloc'd _lhashseed=(double)ihigh; ENDVERBATIM } : unable to get the drand here to recognize the same fseed used in rand FUNCTION vseed () { VERBATIM #ifdef WIN32 if (ifarg(1)) seed=*getarg(1); else { printf("TIME ACCESS NOT PRESENT IN WINDOWS\n"); hxe(); } srand48((unsigned)seed); set_seed(seed); return seed; #else struct timeval tp; struct timezone tzp; if (ifarg(1)) seed=*getarg(1); else { gettimeofday(&tp,&tzp); seed=tp.tv_usec; } srand48((unsigned)seed); set_seed(seed); srandom(seed); ihigh=(unsigned int)seed; return seed; #endif ENDVERBATIM } : mc4seed() sets 'low' bits for mccell_ran4() : opt 2nd arg ihigh is 'high' bits; see (/usr/site/nrniv/nrn/src/ivoc/ivocrand.cpp:85) u_int32_t ihigh_, orig_, ilow_ FUNCTION mc4seed () { VERBATIM double x; u_int32_t idx; unsigned int y; if (!ifarg(1)) { printf("low:%d high:%d\n",(unsigned int)ilow,(unsigned int)ihigh); return (double)ihigh; } else { // u_int32_t idx=(u_int32_t) chkarg(1, 0., 4294967295.); x = *getarg(1); if (x>dmaxuint) {y=(unsigned int)x%4294967296; printf("Warning truncating %20.0f to %d (ilow)\n",x,y);x=y;} ilow = idx = (u_int32_t)x; mcell_ran4_init(idx); } if (ifarg(2)) { x = *getarg(2); if (x>dmaxuint) {y=(unsigned int)x%4294967296; printf("Warning truncating %20.0f to %d (ihigh)\n",x,y);x=y;} ihigh=(u_int32_t) x; } return (double)ihigh; ENDVERBATIM } : from Numerical Recipes in C FUNCTION gammln (xx) { VERBATIM { double x,tmp,ser; static double cof[6]={76.18009173,-86.50532033,24.01409822,-1.231739516,0.120858003e-2,-0.536382e-5}; int j; x=_lxx-1.0; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.0; for (j=0;j<=5;j++) { x += 1.0; ser += cof[j]/x; } return -tmp+log(2.50662827465*ser); } ENDVERBATIM } FUNCTION betai(a,b,x) { VERBATIM { double bt; double gammln(),betacf(); if (_lx < 0.0 || _lx > 1.0) {printf("Bad x in routine BETAI\n"); hxe();} if (_lx == 0.0 || _lx == 1.0) bt=0.0; else bt=exp(gammln(_la+_lb)-gammln(_la)-gammln(_lb)+_la*log(_lx)+_lb*log(1.0-_lx)); if (_lx < (_la+1.0)/(_la+_lb+2.0)) return bt*betacf(_la,_lb,_lx)/_la; else return 1.0-bt*betacf(_lb,_la,1.0-_lx)/_lb; } ENDVERBATIM } VERBATIM #define ITMAX 100 #define EPS 3.0e-7 ENDVERBATIM FUNCTION betacf(a,b,x) { VERBATIM { double qap,qam,qab,em,tem,d; double bz,bm=1.0,bp,bpp; double az=1.0,am=1.0,ap,app,aold; int m; void nrerror(); qab=_la+_lb; qap=_la+1.0; qam=_la-1.0; bz=1.0-qab*_lx/qap; for (m=1;m<=ITMAX;m++) { em=(double) m; tem=em+em; d=em*(_lb-em)*_lx/((qam+tem)*(_la+tem)); ap=az+d*am; bp=bz+d*bm; d = -(_la+em)*(qab+em)*_lx/((qap+tem)*(_la+tem)); app=ap+d*az; bpp=bp+d*bz; aold=az; am=ap/bpp; bm=bp/bpp; az=app/bpp; bz=1.0; if (fabs(az-aold) < (EPS*fabs(az))) return az; } printf("a or b too big, or ITMAX too small in BETACF"); return -1.; } ENDVERBATIM } FUNCTION symval() { VERBATIM { Symbol *sym; sym = hoc_get_symbol(* hoc_pgargstr(1)); // should do type check eg sym->type == VAR return *(hoc_objectdata[sym->u.oboff]._pval); // is this ._pval safe?? } ENDVERBATIM } :* tval(r,N) , computes t-statistic , r == pearson correlation coefficient, N == size of sample FUNCTION tstat() { VERBATIM double r = fabs(*getarg(1)); double N = *getarg(2); if(N < 2) { printf("tstat ERRA: N must be > 2!\n"); return -1; } return r * sqrt(N-2.)/sqrt(1.0-(r*r)); ENDVERBATIM } FUNCTION tdistrib() { VERBATIM double gammln(); double x = *getarg(1); double dof = *getarg(2); double res = (gammln( (dof+1.0) / 2.0 ) / gammln( dof / 2.0 ) ); double pi = 3.14159265358979323846; res *= (1.0 / sqrt( dof * pi ) ); res *= pow((1 + x*x/dof),-1.0*((dof+1.0)/2.0)); return res; ENDVERBATIM } : based on code from http://www.psychstat.missouristate.edu/introbook/rdist.htm : return estimate of critical correlation coefficient (r) that >= to it would be considered significant : at a given sample size of \$1 (N) . probability of null hypothesis p < 0.05 : iff \$2 == 1, probability of null hypothesis p < 0.01 ( r will then be higher ) : null hypothesis == correlation between variables , r , == 0 : this function uses a lookup table so is limited in estimates -- max sample size , N, == 1000 , but : if your sample is bigger and your r value is also bigger, the correlation is significant because : in general, as the sample size , N , increases, the r value above which is considered significant decreases FUNCTION rcrit () { VERBATIM double rtbl[68][3]={//for each row: index0==N, index1==p0.05 value of r, index2==p0.01 value of r 1 , 0.997 , 0.999 , 2 , 0.950 , 0.990 , 3 , 0.878 , 0.959 , 4 , 0.811 , 0.917 , 5 , 0.754 , 0.874 , 6 , 0.707 , 0.834 , 7 , 0.666 , 0.798 , 8 , 0.632 , 0.765 , 9 , 0.602 , 0.735 , 10 , 0.576 , 0.708 , 11 , 0.553 , 0.684 , 12 , 0.532 , 0.661 , 13 , 0.514 , 0.641 , 15 , 0.482 , 0.606 , 16 , 0.468 , 0.590 , 17 , 0.456 , 0.575 , 18 , 0.444 , 0.561 , 19 , 0.433 , 0.549 , 20 , 0.423 , 0.537 , 21 , 0.413 , 0.526 , 22 , 0.404 , 0.515 , 23 , 0.396 , 0.505 , 24 , 0.388 , 0.496 , 25 , 0.331 , 0.487 , 26 , 0.374 , 0.478 , 27 , 0.367 , 0.470 , 28 , 0.361 , 0.463 , 29 , 0.355 , 0.456 , 30 , 0.349 , 0.449 , 31 , 0.344 , 0.442 , 32 , 0.339 , 0.436 , 33 , 0.334 , 0.430 , 34 , 0.329 , 0.424 , 35 , 0.325 , 0.418 , 36 , 0.320 , 0.413 , 37 , 0.316 , 0.408 , 38 , 0.312 , 0.403 , 39 , 0.308 , 0.398 , 40 , 0.304 , 0.393 , 41 , 0.301 , 0.389 , 42 , 0.297 , 0.384 , 43 , 0.294 , 0.380 , 44 , 0.291 , 0.376 , 45 , 0.288 , 0.372 , 46 , 0.284 , 0.368 , 47 , 0.281 , 0.364 , 48 , 0.279 , 0.361 , 58 , 0.254 , 0.330 , 63 , 0.244 , 0.317 , 68 , 0.235 , 0.306 , 73 , 0.227 , 0.296 , 78 , 0.220 , 0.286 , 83 , 0.213 , 0.278 , 88 , 0.207 , 0.270 , 93 , 0.202 , 0.263 , 98 , 0.195 , 0.256 , 123 , 0.170 , 0.230 , 148 , 0.159 , 0.210 , 173 , 0.148 , 0.194 , 198 , 0.138 , 0.181 , 298 , 0.113 , 0.148 , 398 , 0.098 , 0.128 , 498 , 0.088 , 0.115 , 598 , 0.080 , 0.105 , 698 , 0.074 , 0.097 , 798 , 0.070 , 0.091 , 898 , 0.065 , 0.086 , 998 , 0.062 , 0.081 }; double N; int get99, i , tablen, df; N = *getarg(1); get99 = ifarg(2) ? (int)*getarg(2) : 0; tablen = 68; if(N < 3){ printf("rcrit ERRA: N must be >= 3!\n"); return -1.0; } if( N > 998+2 ) { printf("rcrit WARNA: Using N=1000 as estimate.\n"); N = 998+2; } df = (int)N - 2; for(i=0;i df) { if(get99) return rtbl[i-1][2] + ((rtbl[i][2] - rtbl[i-1][2])*((df - rtbl[i-1][0] )/(rtbl[i][0] - rtbl[i-1][0]))); else return rtbl[i-1][1] + ((rtbl[i][1] - rtbl[i-1][1])*((df - rtbl[i-1][0] )/(rtbl[i][0] - rtbl[i-1][0]))); } } return -1.0; ENDVERBATIM }