% sigma_error.m show error for different methods 
function sigma_error
  format compact;
  sprintf('output file sigma_error_m.out being written \n')
  fid = fopen('sigma_error_m.out', 'w');
  n = 100;
  x=zeros(1,n); % faster if declare storage
  fprintf(fid,'sigma_error.m running, n=%d \n', n);

  for i=1:n
    x(i) = 1.0*i;
  end;
  fprintf(fid,'x 1 .. 100 \n');

  sigma1 = stddeva(n, x);
  fprintf(fid,'sum x^2 - (sum x)^2/n  sigma1 = %g \n', sigma1);

  sigma2 = stddevb(n, x);
  fprintf(fid,'sum(x-mean)2/n  sigma2 = %g \n', sigma2);

  for i=1:n
    x(i) = 10.0*i;
  end;
  fprintf(fid,'x 10 .. 1000 \n');

  sigma1 = stddeva(n, x);
  fprintf(fid,'sum x^2 - (sum x)^2/n  sigma1 = %g \n', sigma1);

  sigma2 = stddevb(n, x);
  fprintf(fid,'sum(x-mean)2/n  sigma2 = %g \n', sigma2);

  for i=1:n
    x(i) = 100.0*i;
  end;
  fprintf(fid,'x 100 .. 10000 \n');

  sigma1 = stddeva(n, x);
  fprintf(fid,'sum x^2 - (sum x)^2/n  sigma1 = %g \n', sigma1);

  sigma2 = stddevb(n, x);
  fprintf(fid,'sum(x-mean)2/n  sigma2 = %g \n', sigma2);

  for i=1:n
    x(i) = 1000.0*i;
  end;
  fprintf(fid,'x 1000 .. 100000 \n');

  sigma1 = stddeva(n, x);
  fprintf(fid,'sum x^2 - (sum x)^2/n  sigma1 = %g \n', sigma1);

  sigma2 = stddevb(n, x);
  fprintf(fid,'sum(x-mean)2/n  sigma2 = %g \n', sigma2);

  for i=1:n
    x(i) = 100000000.0*i;
  end;
  fprintf(fid,'x 100000000 .. 10000000000 \n');

  sigma1 = stddeva(n, x);
  fprintf(fid,'sum x^2 - (sum x)^2/n  sigma1 = %g \n', sigma1);

  sigma2 = stddevb(n, x);
  fprintf(fid,'sum(x-mean)2/n  sigma2 = %g \n', sigma2);

  for i=1:n
    x(i) = 100000000.0+i;
  end;
  fprintf(fid,'x 100000001 .. 100000100 \n');

  sigma1 = stddeva(n, x);
  fprintf(fid,'sum x^2 - (sum x)^2/n  sigma1 = %g \n', sigma1);

  sigma2 = stddevb(n, x);
  fprintf(fid,'sum(x-mean)2/n  sigma2 = %g \n', sigma2);
  return; % from sigma_error
  
  function val=stddeva(nn, xx)
    sumx2 = 0.0;
    sumx = 0.0;
    for i=1:nn;
      sumx = sumx + xx(i);
      sumx2 = sumx2 + xx(i)*xx(i);
    end;
    val = sqrt((sumx2 - sumx*sumx/nn)/nn);
    return;
  end % stddeva

  function val=stddevb(nn, xx)
    mean = 0.0;
    sum2 = 0.0;
    for i=1:nn;
      mean = mean + xx(i);
    end;
    mean = mean/nn;
    for i=1:n
      sum2 = sum2 +(xx(i)-mean)*(xx(i)-mean);
    end;
    val = sqrt(sum2/nn);
    return ;
  end % stddevb
end % sigma_error.m