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