Why is my code just filling row 100 values of a matrix?

My code below should fill the matrix q with the value h(when the while loop breaks for each m) however it seems that most of the entries are not filled just left as zero, possibly just the diagonal entries.

I want to run the smaller for loop(for -0.05:0.01:0.05) 100 times and then take the average break time for each gradient.

Would appriciate any help with this.

function [t,I] = Gillespie(beta,gamma,I0,N0,MaxTime)


% Sets up default parameters if necessary.
if nargin == 0
   beta=1/45;
   gamma=1/50;
   mu=0;
   N0=10000;
   I0=ceil(N0*(1-(gamma + mu)/beta));
   MaxTime=12*365;
end


% The main iteration 

[t, po]=Gill_Stochast([0 MaxTime],[I0],[beta gamma N0]);

T=t/365; I=po;

% plots the graphs with scaled colours
subplot(2,1,1)
h=plot(T,I,'-r');
xlabel 'Years';
ylabel 'Infectious'



function [T,P]=Gill_Stochast(Time,Initial,Parameters)


for i = 1:100
a=1;
q = zeros([100,11]);
z = zeros([1,11]);
for m = -0.05:0.01:0.05 %run code for varying gradients of treatment levels

x = [1 2 3 4 5 6 7 8 9 10 11 12];
y = (m*x) + 0.7 - (5.5*m) ; 


subplot(2,1,2)
plot(x,y);

xlabel 'year';
ylabel 'treatment level';
r =365*[1 2 3 4 5 6 7 8 9 10 11 12 13 14];
I=Initial(1);
T=0; P(1)=I;
b=1;
k=1;
while (T(k)<Time(2)&& I>0);

    [step,new]=Iterate(I,Parameters);
    k=k+1;
    T(k)= T(k-1) + step;   %update with time step from iterate
	if     T(k) > r(b) && T(k-1) < r(b) %if time hits 1,2,3 etc years then treat a proportion y(b) of people
           I = I-floor(I*y(b));
           T(k)=r(b);
	       b=b+1;
     
           P(k)=I;
    else
           I = new;
           P(k)=I;
    
    end
    
	
	if I <= 0   %when Infected go to zero then store time as h
        h=T(k);
		break
    else 
        h = 12*365;
	end

    if k>=length(T)  %doubles vector length to skip array list doing so
        T(k*2)=0;
        P(k*2)=0;
    end
end

T=T(1:k); P=P(1:k);

q(i,a)= h;  %updates matrix value (i,a) to be the time that the value I goes to zero for each m (-0.05 corresponds to a=1, -0.04 to a = 2 etc)
z(a)= m;
a = a+1;

end


end

q=q/365; %gives the array q in years rather than days
plot(z,mean(q)); %plots gradients against the mean of each column of q


% Do the actual iteration step
function [step, new_value]=Iterate(I, Parameters)

beta=Parameters(1); gamma=Parameters(2); N=Parameters(3);

RInfect = beta*(N-I)*I/N;
RRecov = gamma*I;

r1=rand(1,1);
r2=rand(1,1);

step = -log(r2)/(RInfect+RRecov);

if r1<RInfect/(RInfect+RRecov)
    new_value=I+1;  % do infection
else
    new_value=I-1;  % do recovery
end

[Geoff Hayes]

reuben - well what is h being set to? If it is

h=T(k);

then what is T(k)? Note that q is 100x11 matrix but you always update it as

q(a) = h;

Is this intentional? Without having access to the inputs (to this function) it will be difficult for us to reproduce your error. Have you tried stepping through the code with the MATLAB debugger?

[Jan]

Start with cleaning up the code. A proper indentation is applied by pressing Ctrl-a Ctrl-i. Remove the commented code, because it is confusing only. Insert meaningful comments, which explain for each line or block, what it is intented. Without any meaningful comments, it is impossible to guess the purpose of the code. In consequence there is no chance to help you.

[reuben crockett]

I have updates some of the things you asked for and included the whole code so is easier to understand. The code runs a gillespie algorithm with increasing or decreasing treatment at each year.

Answer 1

[Geoff Hayes]

reuben - the problem is with these lines of code

for i = 1:100
    a=1;
    q = zeros([100,11]);
    z = zeros([1,11]);

On every iteration of this loop, you are resetting the q matrix to be all of zeros...and so are losing all information from the previous iteration. Instead do

q = zeros([100,11]);
for i = 1:100
    a=1;
    z = zeros([1,11]);

[Jan]

+1. By the way, I'd prefer zeros(100, 11) which avoids the creation of the vector [100,11] only to split it into elements inside the zeros() function afterwards. But this might be a question of taste.