2012-PHYS-339 Statistics / Geiger

From McGill University Physics Department Technical Services Wiki

News

Files

• Geiger handout (http://www.ugrad.physics.mcgill.ca/resources/339/geiger/geiger.pdf) - PDF version of the lab handout
• Geiger Supplement (http://www.ugrad.physics.mcgill.ca/resources/339/geiger/geiger-supplement.pdf) - A bunch of stuff which might be interesting or useful
• sample.data (http://www.ugrad.physics.mcgill.ca/resources/339/geiger/sample.data) - The self-consistent sample data which is characterized in the appendix of the handout
• read_geiger.m (http://www.ugrad.physics.mcgill.ca/resources/339/geiger/read_geiger.m) - MATLAB function to load geiger data
• old-sample.data (http://www.ugrad.physics.mcgill.ca/resources/339/geiger/old-sample.data) - The imfamous sample data as collected by defective acquisition software - can you figure out what went wrong?

Guide to running the provided software

The geiger software is installed on all 339 workstations, and is in the path. There is actually a little bit of trickery going on. The actually geiger software runs on the beige machine, while the command geiger on the workstation is a shell script to SSH in to the beige machine and run the real command in the appropriate subdirectory.

If this is the first time you run the geiger program, you may get asked if you really want to connect ... you do, type "yes".

mark.orchard-webb@ws339-17:~$ geiger
The authenticity of host 'lm339-17 (132.206.252.107)' can't be established.
RSA key fingerprint is 9d:52:72:9d:09:d4:5d:36:a3:f1:b0:f1:e9:67:a9:1f.
Are you sure you want to continue connecting (yes/no)?

If you have not set up public/private keys for SSH between machines, you will get asked for your password each time you run geiger. This is boring. If you have no idea what I'm talking about regarding public/private keys, just run "autossh" once from the command-line ... after this no password should be required.

mark.orchard-webb@ws339-17:~$ autossh 
Generating public/private rsa key pair.
Created directory '/home/morcha/.ssh'.
Your identification has been saved in /home/morcha/.ssh/id_rsa.
Your public key has been saved in /home/morcha/.ssh/id_rsa.pub.
The key fingerprint is:
a9:2b:af:b3:f5:99:32:50:ac:92:dc:15:f5:ee:36:a0 mark.orchard-webb@ws339-17
The key's randomart image is:
+--[ RSA 2048]----+
|       ..        |
|      .  .       |
|     . .  .      |
|      +  o       |
| . o +  S .      |
|  + +  o o       |
|   . .E   +      |
|    o.oo + .     |
|    o*oo+        |
+-----------------+

Responding to all those questions from geiger is also boring. You can run with no questions asked, which is very handy for running from a shell script.

With the Geiger Counter turned on and the gamma source placed close to the detector, try the following:

mark.orchard-webb@artemis:/mnt/data/orchard/geiger$ geiger -a -r 100 -i 100 -p 0.05 -o geiger.data

You should get a window which after a few seconds looks something like this:

Structure of data

Check the documentation in the source tar file

christian.voyer@hades:~$ man ./geiger.5
Reformatting geiger.5, please wait...

Reads

GEIGER(5)                   Lab Programmers Manual                   GEIGER(5)

NAME
       geiger - a data collection tool - the output data format

DESCRIPTION
       The  data  saved by the geiger program consists of two major parts, the
       header and the data-body.

HEADER
       The header is a single line consisting of 4 numbers:  replicas,  inter-
       vals,  columns  and  period.  replicas, intervals and columns are inte-
       gers, period is a floating point number.  The definitions of the param-
       eters are as follows:

       replicas
              The number of replicas requested.

       intervals
              The number of intervals per replica.

       columns
              The number of columns in the histogram.

       period The period in seconds requested by the user.  If the true period
              differed from this because of quantization this a user  problem.

DATA-BODY
       The  data-body  section  consists of replicas rows, each row containing
       columns numbers, separated by tab characters.  The sum of numbers in  a
       given  row  should  equal intervals.  The number in the Nth column of a
       replica represents the frequency of observation of N counts per  period
       during the replica measurement.

SEE ALSO
       geiger(1)

BUGS
       The period problem mentioned above.

AUTHOR
       Mark Orchard-Webb (orchard@physics.mcgill.ca)

Chi-Square, χ², distribution

The following is lifted from Appendix C4 of Bevington.

The probability distribution P_x(x²,ν) for χ² is given by

The probability of observing a value of chi-square larger than χ² for a random sample of N observations with ν degrees of freedom is the integral of this probability P_χ(χ²,ν).

Values of reduced chi-square corresponding to various values of the integral probability P_χ(χ²,ν). of exceeding χ² in a measurement with \nu degrees of freedom are tabulated in Table C-4 for ν ranging from 1 to 200. These values were calculated with the subroutine PCHISQ of Program 10-1. The functional dependence of P_χ(χ²,ν). corresponding to representative values of ν is graphed in Figure C-4 as a smooth variation with the reduced chi-square .

The function which is tabulated and graphed is the shaded area under the tail of the probability curve for values larger than χ² as indicated.

Figure C-4

Table C-4

	P	0.99	0.98	0.95	0.90	0.80	0.70	0.60	0.50	0.40	0.30	0.20	0.10	0.05	0.02	0.01	0.001
ν		======	======	======	======	======	======	======	======	======	======	======	======	======	======	======	======
1		0.00016	0.00063	0.00393	0.0158	0.0642	0.148	0.275	0.455	0.708	1.074	1.642	2.706	3.842	5.412	6.635	10.827
2		0.0100	0.0202	0.0515	0.105	0.223	0.357	0.511	0.693	0.916	1.204	1.609	2.303	2.996	3.912	4.605	6.908
3		0.0383	0.0617	0.117	0.195	0.335	0.474	0.623	0.789	0.982	1.222	1.547	2.084	2.605	3.279	3.782	5.422
4		0.0742	0.107	0.178	0.266	0.412	0.549	0.688	0.839	1.011	1.220	1.497	1.945	2.372	2.917	3.319	4.617
5		0.111	0.151	0.229	0.322	0.469	0.600	0.731	0.870	1.026	1.213	1.458	1.847	2.214	2.678	3.017	4.103
-
6		0.145	0.189	0.273	0.367	0.512	0.638	0.762	0.891	1.035	1.205	1.426	1.774	2.099	2.506	2.802	3.743
7		0.177	0.224	0.309	0.405	0.546	0.667	0.785	0.907	1.040	1.198	1.400	1.717	2.010	2.375	2.639	3.475
8		0.206	0.254	0.342	0.436	0.574	0.691	0.803	0.918	1.044	1.191	1.379	1.670	1.938	2.271	2.511	3.266
9		0.232	0.281	0.370	0.463	0.598	0.710	0.818	0.927	1.046	1.184	1.360	1.632	1.880	2.187	2.407	3.097
10		0.256	0.306	0.394	0.486	0.618	0.727	0.830	0.934	1.047	1.178	1.344	1.599	1.831	2.116	2.321	2.959
-
11		0.277	0.328	0.416	0.507	0.635	0.741	0.840	0.940	1.048	1.173	1.330	1.570	1.789	2.056	2.248	2.842
12		0.298	0.348	0.435	0.525	0.651	0.753	0.849	0.945	1.049	1.168	1.318	1.546	1.752	2.005	2.185	2.742
13		0.317	0.366	0.453	0.542	0.664	0.764	0.856	0.949	1.049	1.163	1.307	1.524	1.720	1.959	2.130	2.655
14		0.333	0.384	0.469	0.556	0.676	0.773	0.863	0.953	1.049	1.159	1.297	1.505	1.692	1.919	2.082	2.580
15		0.349	0.399	0.484	0.570	0.687	0.782	0.869	0.956	1.049	1.155	1.287	1.487	1.666	1.884	2.039	2.462
-
16		0.363	0.413	0.498	0.582	0.697	0.789	0.874	0.959	1.049	1.151	1.279	1.471	1.644	1.852	2.000	2.453
17		0.377	0.427	0.510	0.593	0.706	0.796	0.879	0.961	1.048	1.148	1.272	1.457	1.623	1.823	1.965	2.222
18		0.390	0.439	0.522	0.604	0.714	0.802	0.883	0.963	1.048	1.145	1.264	1.444	1.604	1.797	1.934	2.351
19		0.402	0.451	0.532	0.613	0.722	0.808	0.887	0.965	1.048	1.142	1.258	1.432	1.587	1.772	1.895	2.016
20		0.413	0.462	0.543	0.622	0.729	0.813	0.891	0.967	1.048	1.139	1.252	1.421	1.571	1.751	1.878	2.266
-
22		0.434	0.482	0.561	0.638	0.742	0.823	0.897	0.970	1.047	1.134	1.241	1.401	1.542	1.712	1.831	2.194
24		0.452	0.500	0.577	0.653	0.753	0.831	0.902	0.972	1.046	1.129	1.231	1.383	1.517	1.678	1.791	2.133
26		0.469	0.516	0.592	0.665	0.762	0.838	0.907	0.975	1.045	1.125	1.223	1.368	1.496	1.648	1.756	2.079
28		0.485	0.530	0.605	0.676	0.771	0.845	0.911	0.976	1.045	1.121	1.215	1.354	1.476	1.622	1.724	2.032
30		0.498	0.544	0.617	0.687	0.779	0.850	0.915	0.978	1.044	1.118	1.208	1.342	1.459	1.599	1.697	1.990
-
32		0.511	0.556	0.627	0.696	0.786	0.855	0.918	0.979	1.043	1.115	1.202	1.331	1.444	1.578	1.671	1.953
34		0.523	0.567	0.637	0.705	0.792	0.860	0.921	0.981	1.042	1.112	1.196	1.321	1.429	1.559	1.649	1.919
36		0.534	0.577	0.646	0.712	0.798	0.864	0.924	0.982	1.042	1.109	1.191	1.312	1.417	1.541	1.628	1.889
38		0.545	0.587	0.655	0.720	0.804	0.868	0.926	0.983	1.041	1.107	1.186	1.303	1.405	1.526	1.610	1.861
40		0.554	0.596	0.663	0.726	0.809	0.872	0.928	0.983	1.041	1.104	1.182	1.295	1.394	1.511	1.592	1.835
-
42		0.563	0.604	0.670	0.733	0.813	0.875	0.930	0.984	1.040	1.102	1.178	1.288	1.384	1.498	1.576	1.812
44		0.572	0.612	0.677	0.738	0.818	0.878	0.932	0.985	1.040	1.100	1.174	1.281	1.375	1.485	1.562	1.790
46		0.580	0.620	0.683	0.744	0.822	0.881	0.934	0.986	1.039	1.098	1.170	1.275	1.366	1.473	1.548	1.770
48		0.587	0.627	0.690	0.749	0.826	0.884	0.936	0.986	1.038	1.096	1.167	1.269	1.358	1.463	1.535	1.751
50		0.594	0.633	0.695	0.754	0.829	0.886	0.937	0.987	1.038	1.095	1.163	1.263	1.350	1.452	1.523	1.733
-
60		0.625	0.662	0.720	0.774	0.844	0.897	0.944	0.989	1.036	1.087	1.150	1.240	1.318	1.410	1.473	1.660
70		0.649	0.684	0.739	0.791	0.856	0.905	0.949	0.991	1.034	1.081	1.139	1.222	1.293	1.377	1.435	1.605
80		0.669	0.703	0.755	0.803	0.865	0.912	0.952	0.992	1.032	1.077	1.130	1.207	1.274	1.351	1.404	1.560
90		0.686	0.718	0.768	0.814	0.873	0.917	0.956	0.993	1.031	1.073	1.123	1.195	1.257	1.330	1.379	1.525
100		0.701	0.732	0.779	0.824	0.880	0.921	0.958	0.993	1.030	1.069	1.117	1.185	1.244	1.312	1.358	1.495
-
120		0.724	0.753	0.798	0.839	0.890	0.928	0.962	0.995	1.028	1.064	1.107	1.169	1.221	1.283	1.325	1.447
140		0.743	0.770	0.812	0.850	0.898	0.934	0.965	0.995	1.026	1.059	1.099	1.156	1.204	1.261	1.299	1.410
160		0.759	0.784	0.824	0.860	0.905	0.939	0.968	0.996	1.024	1.056	1.093	1.146	1.191	1.243	1.278	1.381
180		0.771	0.796	0.833	0.868	0.910	0.942	0.970	0.996	1.023	1.052	1.088	1.137	1.179	1.228	1.262	1.358
200		0.782	0.806	0.841	0.874	0.915	0.945	0.972	0.997	1.022	1.050	1.083	1.130	1.170	1.216	1.247	1.338

Various codes

The fortran code for PCHISQ is shown below

      FUNCTION PCHISQ (CHISQ, NFREE)
      DOUBLE PRECISION Z, TERM, SUM
 11   IF (NFREE) 12, 12, 14
 12   PCHISQ = 0.
      GO TO 60
 14   FREE = NFREE
      Z = CHISQ*FREE/2
      NEVEN = 2*(NFREE/2)
      IF (NFREE - NEVEN) 21, 21, 41
C
C     NUMBER OF DEGREES OF FREEDOM IS EVEN
C     
 21   IMAX = NFREE/2
      TERM = 1.
      SUM = 0.
 31   DO 34 I=1, IMAX
      FI = I
      SUM = SUM + TERM
 34   TERM = -TERM * Z/FI
 35   PCHISQ = SUM * DEXP(-Z)
      GOTO 60
C     
C     NUMBER OF DEGREES OF FREEDOM IS ODD
C     
 41   IF (Z - 25) 44, 44, 42
 42   Z = CHISQ * (FREE-1.)/2.
      GO TO 21
 44   PWR = FREE/2
      TERM = 1.
      SUM = TERM/PWR
 51   DO 56 I=1, 1000
      FI = I
      TERM = TERM * Z/FI
      SUM = SUM + TERM/(PWR+FI)
 55   IF (DABS(TERM/SUM) - 0.00001) 57, 57, 56
 56   CONTINUE
 57   PCHISQ = 1. - (Z**PWR)*SUM/GAMMA(PWR)
 60   RETURN
      END

This is a little too ugly for me, so I rewrote it in C as follows

double pchisq (double chisq, int nfree)
{
  double free, z, pwr, term, sum, fi;
  int i, imax;
  if (nfree < 1) return 0;
  free = nfree;
  z = chisq * free / 2;
  if (nfree % 2) {		/* nfree is odd */
    if (z > 25) {
      z = chisq * (free-1.)/2.;
    } else {
      pwr = free/2;
      term = 1.;
      sum = term/pwr;
      for (i = 0; i<1000;) {
	i++;
	fi = i;
	term = -term * z/fi;
	sum = sum + term/(pwr+fi);
	if (fabs(term/sum) <= 0.00001) break;
      }
      return 1. - pow(z,pwr)*sum/exp(lgamma(pwr));
    }
  }
  imax = nfree/2;
  term = 1.;
  sum = 0.;
  for (i = 0; i < imax;) {
    i++;
    fi = i;
    sum = sum + term;
    term = term * z/fi;
  }
  return sum * exp(-z);
}

Using the following main

int main (int argc, char **argv)
{
  double rchisq;
  int nu;
  sscanf (argv[1], "%d", &nu);
  for (rchisq = 0; rchisq < 1.7; rchisq+=0.01) {
    fprintf (stdout, "%g %g\n", rchisq, pchisq (rchisq, nu));
  }
  return 0;
}

And the following bit of shell script

#! /bin/sh

# prepare the graph
cat << EOF > pchisq.gnu
set xlabel '{/Symbol C^2/n}'
set ylabel 'P_{/Symbol C}({/Symbol C}^2,{/Symbol n})'
EOF

for nu in 1 2 3 5 10 20 30
do
    ./pchisq $nu > pchisq-$nu.data
    if test $nu = 1
    then 
      echo -n "plot " >> pchisq.gnu
    else
      echo -n ", " >> pchisq.gnu
    fi
    echo -n "'pchisq-$nu.data' w l title '{/Symbol n} = $nu'" >> pchisq.gnu
done

I am able to recreate Bevingtons figure C-4.

Table C-4 was created with the following code, an then hand edited for the variable precision in the upper left corner. As can be seen it started out as an elegant binary search and degenerated in to a brute force linear search as time ran out.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

double pchisq (double chisq, int nfree)
{
  double free, z, pwr, term, sum, fi;
  int i, imax;
  if (nfree < 1) return 0;
  free = nfree;
  z = chisq * free / 2;
  if (nfree % 2) {		/* nfree is odd */
    if (z > 25) {
      z = chisq * (free-1.)/2.;
    } else {
      pwr = free/2;
      term = 1.;
      sum = term/pwr;
      for (i = 0; i<5000;) {
	i++;
	fi = i;
	term = -term * z/fi;
	sum = sum + term/(pwr+fi);
	if (fabs(term/sum) <= 0.00001) break;
      }
      return 1. - pow(z,pwr)*sum/exp(lgamma(pwr));
    }
  }
  imax = nfree/2;
  term = 1.;
  sum = 0.;
  for (i = 0; i < imax;) {
    i++;
    fi = i;
    sum = sum + term;
    term = term * z/fi;
  }
  return sum * exp(-z);
}

double binsearch(int nu, double target, double low, double high)
{
  static int depth = 1;
  double mid = 0.5 * (low + high);
  double P = pchisq(mid, nu);
  depth++;
  if (depth > 20000) {
    fprintf (stderr, "Tooo deeeeeep. nu = %d, target = %f, low = %f, high = %f, P = %f Abort\n", nu, target, low, high, P);
    return -1;
  }
  if ((fabs(target - P)/target) < 0.0001) return mid;
  if (P > target) return binsearch(nu, target, mid, high);
  if (P < target) return binsearch(nu, target, low, mid);
  fprintf (stderr, "No can not get here!!!!!\n");
  exit(1);
}

const char *pretty(double x)
{
  static char buf[16];
  int scale = log(x) / log(10);
  sprintf (buf, "%5.3f", x);
  return buf;
}

double ugly (int nu, double target, double start)
{
  int i;
  for (i = 10000*start; ; i++) {
    if (pchisq(1e-4*i, nu) < target) return 1e-4*i;
  }
  fprintf (stderr, "No can not get here!!!!!\n");
}

int main (int argc, char **argv)
{
  char *Ps[] = {"0.99", "0.98","0.95", "0.90", "0.80","0.70","0.60","0.50","0.40","0.30","0.20","0.10","0.05","0.02","0.01","0.001", 0};
  double *P;
  int ncols, i, nu;
  fprintf (stdout, "<table><tr><td></td><td>P</td>");
  for (i = 0; Ps[i]; i++) fprintf (stdout, "<td>%s</td>", Ps[i]);
  fprintf (stdout, "</tr>\n");
  fprintf (stdout, "<tr><td><math>\\nu</math></td><td></td>");
  ncols = i;			/* take advantage of knowing how many columns */
  P = malloc (ncols*sizeof(double));
  for (i = 0; i < ncols; i++) {
    sscanf (Ps[i], "%lf", &P[i]); /* string to double */
    fprintf (stdout, "<td>======</td>");
  }
  fprintf (stdout, "</tr>\n");  
  for (nu = 1; nu < 10; nu++) {
    fprintf (stdout, "<tr><td>%d</td><td></td>", nu);
    for (i = 0; i < ncols; i++) {
      fprintf (stdout, "<td>%s</td>", pretty(binsearch(nu,P[i],0,12/sqrt(nu))));
    }
    fprintf (stdout, "</tr>\n");
  }
  for (nu = 10; nu < 15; nu++) {
    fprintf (stdout, "<tr><td>%d</td><td></td>", nu);
    for (i = 0; i < ncols; i++) {
      fprintf (stdout, "<td>%s</td>", pretty(binsearch(nu,P[i],0.25,10/sqrt(nu))));
    }
    fprintf (stdout, "</tr>\n");
  }
  for (nu = 15; nu < 20; nu++) {
    fprintf (stdout, "<tr><td>%d</td><td></td>", nu);
    for (i = 0; i < ncols; i++) {
      fprintf (stdout, "<td>%s</td>", pretty(ugly(nu,P[i],0.2)));
    }
    fprintf (stdout, "</tr>\n");
  }
  for (nu = 20; nu < 50; nu+=2) {
    fprintf (stdout, "<tr><td>%d</td><td></td>", nu);
    for (i = 0; i < ncols; i++) {
      fprintf (stdout, "<td>%s</td>", pretty(ugly(nu,P[i],0.4)));
    }
    fprintf (stdout, "</tr>\n");
  }
  for (nu = 50; nu < 100; nu+=10) {
    fprintf (stdout, "<tr><td>%d</td><td></td>", nu);
    for (i = 0; i < ncols; i++) {
      fprintf (stdout, "<td>%s</td>", pretty(ugly(nu,P[i],0.59)));
    }
    fprintf (stdout, "</tr>\n");
  }
  for (nu = 100; nu < 201; nu+=20) {
    fprintf (stdout, "<tr><td>%d</td><td></td>", nu);
    for (i = 0; i < ncols; i++) {
      fprintf (stdout, "<td>%s</td>", pretty(ugly(nu,P[i],0.7)));
    }
    fprintf (stdout, "</tr>\n");
  }
  fprintf (stdout, "</table>\n");
  return 0;
}

Supplemental toys

• chisq.c (http://www.ugrad.physics.mcgill.ca/resources/339/geiger/chisq.c) - This is an implementation of table C4. in Bevington. There are two functions of interest in this file:

double chisq_five (int dof)

returns the reduced Χ² value which only 5% of replicas having dof are expected to exceed.

double chisq_ten (int dof)

returns the reduced Χ² value which only 10% of replicas having dof are expected to exceed.

These functions are written in such a way that they can easily be pasted into another piece of code.

Code to generate Gnuplot script

Here is some delicious copy-pasta which will generate a gnuplot script, which plots nice smooth poisson and gaussian functions. It should look something like

The function expects 7 arguments:

• graphname: this is a string, it will be the name of the gnuplot script to be produced.
• dataname: this is the name of the data file to be plotted, it should contain the usual 3 columns.
• A_p: The normalization constant for the poisson curve (intervals)
• mu_p: The mean of the poisson curve (replica mean)
• A_g: The normalization constant of the gaussian curve (GC * intervals, where GC is great than or equal to 1, due to lost area of gaussian for x < 0). "GC" can be calculated by symbolically evaluating the gaussian integral from 0 to .
• mu_g: The mean of the gaussian curve (replica mean)
• sigma_g: The standard deviation of the gaussian curve (sqrt(replica variance))

int make_gnuplot (const char *graphname, const char *dataname, double A_p, double mu_p, double A_g, double mu_g, double sigma_g)
{
  FILE *graph;

  graph = fopen (graphname, "w"); /* try to open requested filename */
  if (graph == NULL) {		/* if it didn't work ... */
    fprintf (stderr, "Can't open '%s' for write", graphname); /* warn the user */
    exit(1);			/* crash and burn ... "#include <stdlib.h>" to avoid gcc warning */
  }

  /* Dump out the variables */
  fprintf (graph, "A_p = %g\n", A_p);
  fprintf (graph, "mu_p = %g\n", mu_p);
  fprintf (graph, "A_g = %g\n", A_g);
  fprintf (graph, "mu_g = %g\n", mu_g);
  fprintf (graph, "sigma_g = %g\n", sigma_g);

  /* Define the functions */
  fprintf (graph, "poisson(x,A,mu) = A*exp(x*log(mu)-lgamma(x+1)-mu)\n");
  fprintf (graph, "gaussian(x,A,mu,sigma) = A*exp(-0.5*((x-mu)/sigma)**2)/(sqrt(2*pi)*sigma)\n");

  /* Make some pretty labels */
  fprintf (graph, "set xlabel 'Events per 0.1 s interval'\n");
  fprintf (graph, "set ylabel 'Frequency of observation'\n");

  /* The big plot command */
  fprintf (graph, "plot [-1:] '%s' using 1:2:3 title 'sample replica' with error, \\\n", dataname);
  fprintf (graph, "(x<0)?1/0:poisson(x, A_p, mu_p) title 'Poisson' with line, \\\n"); /* See comments on wiki */
  fprintf (graph, "(x<0)?1/0:gaussian(x, A_g, mu_g, sigma_g) title 'Gaussian' with line\n"); /* regarding "(x<0)?1/0:" */
  fclose (graph);
  return 0;
}

Note: portion of code x<0)?1/0:poisson(x, A_p, mu_p) is a dirty (elegant) trick to avoid plotting the curve for domain values less than zero. Basically this reads "plot 1/0 if x < 0, otherwise plot poisson (x, A_p, mu_p)." 1/0 is of course undefined, so gnuplot plots nothing.

Writing out data from a replica to a file

This is a related problem. The following copy-pasta function accepts 4 arguments,

• dataname: this is a string, which is the name of the datafile to be created.
• columns: the number of columns in the replica
• replica: this is the array of frequencies. This bit is tricky! Let's assume you allocated storage for the frequencies in an array as follows:

int data[rows][cols];

You can then pass the n^th replica to this function as

make_graphdata ("sample-replica.data", cols, data[n], colrep);

• colvar: this is the array of column variances to be used as error-bar.

int make_graphdata (const char *dataname, int columns, int *replica, double *colvar)
{
  FILE *datafile;
  int i;

  datafile = fopen (dataname, "w"); /* create a data stream */
  if (datafile == NULL)	{	/* if it failed */
    fprintf (stderr, "Can't open '%s' for write\n", dataname);
    exit(1);			/* going down with the ship  ... "#include <stdlib.h>" to avoid gcc warning */
  }
  for (i = 0; i < columns; i++) {
    if (colvar[i] > 0) {
      fprintf (datafile, "%d\t%d\t%g\n",
	       i, replica[i],
	       sqrt(colvar[i]));
    }
  }
  return 0;
}

Historical Archives

• 2006-PHYS-339 Statistics / Geiger
• 2007-PHYS-339 Statistics / Geiger
• 2008-PHYS-339 Statistics / Geiger
• 2009-PHYS-339 Statistics / Geiger
• 2011-PHYS-339 Statistics / Geiger

Log