1D ising model - Helix-coil transistion

[Päällikkö]

So medres in my code collects several data points at each temperature, say 10000. res, on the other hand, stores an average of these 10000 values to represent the result at each particular temperature. What you want is not only the mean of medres (which is res), but also the standard deviation, sigma (which is an error bar estimate, although quite often people use 2*sigma as their error estimate, or at LHC/CERN, 5*sigma when they announced they had found the Higgs boson). From what I understand, you are trying to achieve a universal way of going about this with your statistics class, which is a good idea and ought to clarify the code.

If you paste your main, or better yet, a minimal example of how you want to use your statistics class (that is, one without interaction with the Particle class etc.), things would probably be easier to debug (if there sitll are bug). This might also clarify to yourself what the limits of your clss design might be (for example, it is not immediately clear to me how you suppose to use the class for multiple temperatures, maybe it is just a dropin for medres, rather than res?).

Finally, I'd recommend using

statistics::statistics() : n(0), sum(0.0), sumsq(0.0) {}

Here it doesn't really matter, but suppose n were const int rather than an int. Only the command I listed would then work (you can initialize with a value, but cannot assign one; same is true if you tried const int x; x = 1; versus const int x = 1; The latter is cheaper too when not talking about consts).

 

[Matt atkinson]

well for example i want to calculate $\sigma$ for each mean value of medres, for each temperature. So inside the loop for temperature, i could put something like this ;

A=medres.getAverage();
B=medres.getSqAverage();
Cnum=medres.getNumber();
Sigma=std::sqrt((B-A*A)/(Cnum));

Not sure if I am understanding it right.

 

[Päällikkö]

That'd be correct, yes, if you define "statistics medres;" and rather than "medres.push_back(particle.mean());" you use "medres.add(particle.mean());"
Sigma should be an array like res (res stores the means, Sigma should store the sigmas at each temperature), or rather res could be replaced by A should A also be made into an array. Given that you compute A, B, C using the methods of statistics, you'd do well to make Sigma a method of statistics, too (say, getError(), or whatever you consider appropriate).

 

[Matt atkinson]

Oh okay thankyou! ill try that when i get back to my code.
Define sigma, in the statistics class then have something like;

sigma=std::sqrt((sumsq-sum*sum)/(n))

i believe this should work
also with

class statistics{
      private:
              int n;
              double sum;
              double sumsq;
              double sigma;
              double medres;
      public:
             statistics();
             int getNumber() const;
             double getAverage() const;
             double getSqAverage() const;
             void add(double x);
};

If i have understood correctly

 

[Matt atkinson]

thankyou for all your help, just checked and all the errors i get are the exact same (~0.00042) have i done something incorrectly?
I think its my definition of medres2.

int main() {
  std::uniform_real_distribution<float> dist{0.f, 1.f};
  std::mt19937                          gen;

  std::vector<float> Ts;
  std::ofstream ofile;
  std::array<float, N>   A,B,C;
  statistics medres2;
  float error;
  ofile.open ("helixfrac.txt");
  for(float T = 420.f; T < 434.f; T+=.025f)
    Ts.push_back(T);
  std::vector<float> res;
  for(auto T : Ts) {
    Particle particle(T);
    std::vector<float> medres;
    medres.reserve(10000);
    for(int i = 0; i < 100000; ++i) {
      particle.choose();
      auto Eb = particle.localEnergy();
      particle.perturb();
      auto Ef = particle.localEnergy();
      if(std::exp(-(Ef-Eb)/T) < dist(gen))
        particle.perturb();
      if(!(i%10) && i > 1000){
        medres.push_back(particle.mean());
        medres2.add(particle.mean());
                            }
    }
    float xx = 0.f;
    for(const auto x : medres)
      xx += x;
    res.push_back(xx/medres.size());
  }
  for(unsigned int i = 0; i < Ts.size(); ++i){
    A[i]=medres2.getSqAverage();
    B[i]=medres2.getAverage();
    C[i]=1/(double) medres2.getNumber();
    std::cout << Ts[i] << "," << (1+res[i])*.5 << std::endl;
    ofile << Ts[i] << "," << (1+res[i])*.5 << "," << std::sqrt((A[i]*C[i]-B[i]*B[i]*C[i])) << std::endl;
    }
    ofile.close();
  std::cout << std::endl;
}

 

[Päällikkö]

It sounds more or less correct. There are caveats in the way that you compute the error (it gives an underestimate), but it is not really your fault, it's a simplification made in your course materials. The reason is that the measurements are highly correlated with each other (the length of the polymer is, perhaps, N = 512, while the !(i%10) means that you are storing an average after every 10 attempts at flipping: Obviously even if all the flips were successful, you'd only move the average computer by Particle::mean() a few percentages at best). If you are interested in better methods, look up block-averaging and other autocorrelation based calculations (the idea is that you wait for as long as there is no correlation in the mean data to take averages). Vary N and the number of flips per temperature (i < 10000) and see how they affect the error estimates.

Finally, there's no reason for you to keep the medres's in your code. I cleaned it up a bit (still not pretty, but now it doesn't have anything that's not needed). I also use a couple of different loop styles, so you can compare and see what you feel comfortable with. For example, stats is an array (std::vector) of Statistics objects, so that each temperature gets its own Statistics. Now I initialize stat as the iterator (if you're not familiar with iterators, just think of it as a pointer) to the first Statistics object in the array, and then increment the iterator at each temperature. When I print out, I use a more traditional loop, where everything is indexed (note, though, that many data types, such as std::list, do not permit direct indexing, i.e. random access, but rather you must traverse through them).

int main() {
  auto dist = std::uniform_real_distribution<float>{0.f, 1.f};
  auto gen  = std::mt19937{};

  auto Ts = std::vector<float>{};
  for(float T = 426.f; T < 428.f; T+=.05f)
    Ts.push_back(T);

  std::vector<Statistics> stats(Ts.size());

  auto stat = stats.begin();
  for(auto T : Ts) {
    auto particle = Particle(T);
    for(int i = 0; i < 10000; ++i) {
      particle.choose();
      auto Eb = particle.localEnergy(); 
      particle.perturb();
      auto Ef = particle.localEnergy();
      if(std::exp(-(Ef-Eb)/T) < dist(gen))
        particle.perturb();
      if(!(i%100) && i > 5000)
        stat->add(particle.mean());
    }
    ++stat;
  }
  for(unsigned int i = 0; i < Ts.size(); ++i) {
    auto mean  = stats[i].getAverage();
    auto stdev = stats[i].getStdDeviation();
    std::cout << Ts[i] << " " 
              << (1+mean)*.5 << " " 
              << (1+mean-2*stdev)*.5 << " " 
              << (1+mean+2*stdev)*.5 << std::endl;
  }
  std::cout << std::endl;
}

You'll want to add the file outstream (i.e. std:fstream) instead of std::cout to write directly to a file. The reason I'm using 2*stdev is that I consider this a better estimate of the error than stdev (and this is indeed more commonly used). You might want to change this.

Finally, as you may notice, I made some modifications to your statistics class:

class Statistics {
public:
  int    getNumber() const;
  float  getAverage() const;
  float  getStdDeviation() const;
  void   add(double x);
private:
  int    n     = 0;
  double sum   = 0.f;
  double sumsq = 0.f;
};

int Statistics::getNumber() const {
  return n;
}

float Statistics::getAverage() const {
  if(n==0) 
    return 0.0;
  return sum/n;
}

float Statistics::getStdDeviation() const {
  if(n==0) 
    return 0.0;
  return std::sqrt(sumsq - sum*sum/n)/n;
}

void Statistics::add(double x) {
   n++;
   sum += x;
   sumsq += x*x;
}
[[SIZE=4]/code][/SIZE]

 

[Matt atkinson]

Wow thank you, your help has been greatly appreciated with this. I learned quite a lot about coding and physics too, while also doing the problem.
A lot of the code you used, i didn't understand at first but your explanations were great and anything else I looked up.

Thanks again!