00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 float mdp_jackboot_plain(float *x, void *a);
00014
00033 class mdp_jackboot {
00034 private:
00035 friend float mdp_jackboot_plain(float *x, void *a) {
00036 int *i_ptr=(int*) a;
00037 return x[*i_ptr];
00038 }
00039 int mdp_jackboot_plain_int;
00040 public:
00041 float (*f)(float*, void *);
00042 int nconf,narg,conf;
00043 float *m;
00044 void* handle;
00045 mdp_jackboot() {
00046 m=0;
00047 dimension(0,0);
00048 }
00050 mdp_jackboot(int nconf_, int narg_=1) {
00051 m=0;
00052 dimension(nconf_, narg_);
00053 }
00054 void dimension(int nconf_, int narg_=1) {
00055 nconf=nconf_;
00056 narg=narg_;
00057 conf=0;
00058 handle=0;
00059 if(m!=0) {
00060 delete[] m;
00061 m=0;
00062 }
00063 if(m==0) m=new float[nconf*narg];
00064 int i,j;
00065 for(i=0; i<nconf; i++)
00066 for(j=0; j<narg; j++)
00067 m[i*narg+j]=0;
00068 }
00069 virtual ~mdp_jackboot() {
00070 delete[] m;
00071 }
00072 float* address(int conf) {
00073 return m+conf*narg;
00074 }
00075 float &operator() (int present_conf, int arg) {
00076 conf=present_conf;
00077 return m[conf*narg+arg];
00078 }
00079 float &operator() (int arg) {
00080 return m[conf*narg+arg];
00081 }
00082 void plain(int i) {
00083 if((i<0) || (i>=narg)) error("incorrect mdp_jackboot.plain() argument");
00084 mdp_jackboot_plain_int=i;
00085 handle=(void*) &mdp_jackboot_plain_int;
00086 f=mdp_jackboot_plain;
00087 }
00088 void makesample(int *p, int nboot) {
00089 int boot, j;
00090 for(boot=0; boot<nboot; boot++)
00091 for(j=0; j<=conf; j++)
00092 p[j+(conf+1)*boot]=(int) ((conf+1)*mdp_random.plain());
00093 }
00094
00095 float mean() {
00096 int i,j;
00097 float tmp;
00098 float *average= new float[narg];
00099 if(conf==0) return (*f)(&m[0],handle);
00100 for(i=0; i<narg; i++) {
00101 average[i]=0;
00102 for(j=0; j<=conf; j++)
00103 average[i]+=m[j*narg+i]/(conf+1);
00104 }
00105 tmp=(*f)(average,handle);
00106 delete[] average;
00107 return tmp;
00108 }
00109 float j_err() {
00110 int i,j,k;
00111 float mean0=this->mean();
00112 float *average= new float[(conf+1)*narg];
00113 float stddev=0;
00114 if(conf<2) return 0;
00115 for(i=0; i<narg; i++)
00116 for(j=0; j<=conf; j++) {
00117 average[j*narg+i]=0;
00118 for(k=0; k<=conf; k++)
00119
00120 if(j!=k) average[j*narg+i]+=(m[k*narg+i]/(nconf-1));
00121 }
00122 stddev=0;
00123 for(j=0; j<=conf; j++)
00124 stddev+=(float) pow((double) (*f)(&(average[j*narg]),handle)-mean0,2.0);
00125 delete[] average;
00126 return (float) sqrt(stddev*conf/(conf+1));
00127 }
00128 float b_err(int nboot=100) {
00129 int i,j,boot,x,y;
00130 float vx, vy, tmp;
00131 float *average=new float[nboot*narg];
00132 int *p=new int[(conf+1)*nboot];
00133 makesample(p,nboot);
00134 if(conf==0) return 0;
00135 for(i=0; i<narg; i++)
00136 for(boot=0; boot<nboot; boot++) {
00137 average[boot*narg+i]=0;
00138 for(j=0; j<=conf; j++)
00139 average[boot*narg+i]+=m[p[j+(conf+1)*boot]*narg+i];
00140 average[boot*narg+i]/=conf+1;
00141 }
00142 for(x=1; x<nboot; x++) {
00143 vx=(*f)(&(average[x*narg]),handle);
00144 for(y=x-1; y>=0; y--) {
00145 vy=(*f)(&(average[y*narg]),handle);
00146 if(vy>vx)
00147 for(i=0; i<narg; i++) {
00148 tmp=average[(y+1)*narg+i];
00149 average[(y+1)*narg+i]=average[y*narg+i];
00150 average[y*narg+i]=tmp;
00151 } else y=-1;
00152 }
00153 }
00154 vx=(*f)(&(average[((int)((float) nboot*0.16))*narg]),handle);
00155 vy=(*f)(&(average[((int)((float) nboot*0.84))*narg]),handle);
00156 delete[] average;
00157 delete[] p;
00158 return (float) (vy-vx)/2.0;
00159 }
00160 };
