How do I loop a fucntion over matrix columns and store the results in a column?

Hi,

I have an issue writing a while loop correctly. I have a variable (R) with values that are stored in a matrix full of columns and a function (dfaedit), which I would like loop over each of those columns. Each of these loops should produce a result (H) which I would like to store in a seperate matrix column as a list of numerical values.

I was wondering if someone could help me write this while-loop correctly?

Thanks!

Answer 1

[Torsten]

If R is your matrix:

[m,n] = size(R);
H = zeros(n,1);
for i = 1:n
  r = R(:,i);
  H(i) = dfaedit(r );
end

[lil brain]

Hi Torsten,

thanks for the quick help. I seem to get an output variable but the values dont seem right to me. The code also seems to throw the following error:

Error using dfaedit_2
Too many input arguments.

Here is the function I am using.

function  [H]=dfaedit(file_name,plot_flag, outfile_flag, out_command_flag)

%  DFA  Calculates and plots detrended fluctuations
%
%   For a 1-d time series in a one column text file named 'file_name',
%	DFA performs a detrended fluctuation analysis on this data. It first integrates
%   the data and then partitions it into different size subsets whose size  
%   varies from MinBox (default = 4) to MaxBox (default = 1/4 the data length). 
%   For the subsets of each size, the best fitting trend lines are found 
%   and their the root mean square residual (Q) is calculated.   The log10
%   of this fluctuation quantity (Q) is plotted against log10 of the points in a subset.
%   The slope of this regression estimates the Hurst exponent (H).  The fractal 
%   dimension can be calculated as 2-H.
%
%   If plot_flag = 1 then a regression plot is created.  If outfile_flag = 1, the file name, H
%   and r^2 are appended to a text file 'dfa_results.txt'.  If out_command_flag = 1, the data 
%   that are the basis for the regression plot are output to the command window.
%
%   Example:  dfa('noise-f.dat', 1, 1, 1)
%
%   I used the Peng's C code algorithm which is available at http://reylab.bidmc.harvard.edu.   
%   It is discussed in Peng et al. Quantification of scaling exponents and crossover  
%   phenomena in  nonstationaryheartbeat time series. Chaos 1995;5:82-87.  I vectorized the 
%   preliminary calculations for MatLab.  But the main loop I did not touch. Still not sure 
%   of the algorithm they are using to calculate the root mean square residual. 
%
%   Syntax:
%   dfa(file_name,plot_flag, outfile_flag, out_command_flag)
%
%   R. C. Schmidt 11-27-01
%-------------------------------------------------------------------------------------
format short g
x = file_name;
numberpoints = length(x);

MinBox = 4;                 % minbox #of points in a box
MaxBox  = .25*numberpoints; % maxbox set 1/4 the data length
BoxSizeDensity = 4;
LogScaleFactor = power(2.0, 1.0/BoxSizeDensity);
index = 1:numberpoints; 
index = reshape(index, length(index), 1);  %change row to column vector

% Preliminary calculations
Sum = cumsum(x);
SumOfSumsSquared = cumsum(Sum.*Sum);
SumOfSums = cumsum(Sum);
SSc = cumsum(index.*Sum);

% now find best fit lines and find fluctuation about the line 
% loop for each box size, from MinBox to MaxBox
iteration = 1;
BoxSize = MinBox;
TempBoxSize = MinBox;
while BoxSize <= MaxBox
	  s = 0.0;
	  r1 = 1./(BoxSize + 1.0);
	  Det = 12.0*r1*r1*r1/(1.0 - r1*r1);
	  vv = BoxSize*(2.0*BoxSize + 1.0)/6.0;

	  for j = 2:numberpoints - BoxSize
	    dhh = SumOfSumsSquared(j + BoxSize) - SumOfSumsSquared(j-1);
	    dhu = SumOfSums(j + BoxSize) - SumOfSums(j-1);
	    dhv = SSc(j + BoxSize) - SSc(j-1) - dhu*j;

	    s = s + dhh - (dhv*dhv +dhu*dhu*vv - dhv*dhu*BoxSize)*Det;
		j = j + 1;
	   end
	
	  den = (numberpoints - BoxSize)*(BoxSize + 1.0);
	  Fluctuation = sqrt(s/den);
	  
      log_points_in_box(iteration,1) = log10(BoxSize);  
	  log_Q(iteration,1) = log10(Fluctuation);
	  show_results(iteration,:) = [ iteration  BoxSize  Fluctuation  log_points_in_box(iteration,1)  log_Q(iteration,1) ];
	  iteration = iteration + 1;

% update the box size 
	  TempBoxSize = TempBoxSize*LogScaleFactor;
	  while round(TempBoxSize) <  BoxSize + 1.0
        TempBoxSize = TempBoxSize*LogScaleFactor;
      end
	  BoxSize = round(TempBoxSize);
end

% got all boxes; now calculate H via trendline
r_trend = corrcoef(log_points_in_box, log_Q);
coefs = polyfit(log_points_in_box, log_Q,1);
r_line = polyval(coefs,log_points_in_box);

% calculate dimension and Hurst
H = coefs(1);
D = 2-H;

%display results in command window
if out_command_flag == 1
disp('         iteration    points       Q       Log10(pnts)   Log10(Q)')
disp(show_results)
end

% display results in a figure
if plot_flag == 1
	figure
	h = axes('Position', [0 0 1 1], 'Visible', 'off');
	axes('Position',[.1 .2 .75 .75])
	plot(log_points_in_box, log_Q,'b+')
	hold on
	plot(log_points_in_box,r_line, 'k')
	ylabel('log10(Q)');
	xlabel('log10(Points in Subset)');
	title(file_name);
	set(gcf, 'CurrentAxes', h);
	str(1) = {[sprintf('%2.5f',coefs(1)),'x + ', sprintf('%2.5f',coefs(2)), ' = y,   r^2 = ', sprintf('%2.2f',r_trend(1,2)^2), ',  n = ', sprintf('%d',length(log_Q)) ]};
	text(.1, .02, str, 'FontSize', 10, 'Color', 'r');
	str(1) = {['H = ', sprintf('%2.3f',H) '  D = ', sprintf('%2.3f',D) ]};
	text(.5, .02, str, 'FontSize', 10, 'Color', 'r');
end

% append results to a file rd_results.txt
if outfile_flag == 1
	fid = fopen('dfa_results.txt','a');
	fprintf(fid,'%s  %4.4f   %2.4f\n', file_name, H, r_trend(1,2)^2);
	status = fclose(fid);
end

Any idea what is going wrong?

Also, I have attached the matrix R in case that helps.

Thanks!

[Torsten]

And you call this function as

[m,n] = size(R);
H = zeros(n,1);
for i = 1:n
    r = R(:,i);
    H(i) = dfaedit(r,1,1,1);
end

?

And what is dfaedit_2 ?

[lil brain]

This works! Very cool and thank you for the support. The dfaedit_2 was just a different version of the same function I used :)