@Markiyan Hirnyk 

Dear Markiyan Hirnyk,

thanks for the answer, and the smart, compact solution proposed. I didn't know about the DirectSearch package.

jon

@jonlg 

Found another typo, but the key issue is not solved yet: Mapleprimes_LQ_model_simulation.mw

Thanks,

jon

@jonlg 

Hi there,

there was still one more typo. This is the new worksheet: Mapleprimes_LQ_model_simulation.mw

I'm still about the correct approach to define the functions, taking into account that I have one varying parameter for the single dose case, and tow of the for the fractionated case.

Thanks,

jon

@Carl Love 

Dear Carl,

thanks for your answer. Sorry for having it in 2D Input :-S I was recommeneded it.

Yes, not closing the quotes was a quite a big error.

On the other hand, it looks like the functions of interest (SHs, STs, SHf, STf) do not get evaluated or I'm lacking to specify t at some point to have them evaluated.

I've also have some concerns about the commands

Ds:=t->tDs:

map(Gf, (D1f^2+D2f^2+2*D1f*D2f*exp(-lambda*T))/tDf^2, T)

They may not be correctly stated. The truth is that for the single dose case 8the first one above, Ds is constant), but I guess the above statement says that it is indeed constant for all t.

This is the new worksheet having corrected the quotes: Mapleprimes_LQ_model_simulation.mw

thanks,

jon

Hi there,

any guess?

Thanks,

jon

@Carl Love 

Dear Carl,

thanks again. Again, you pointed in the right direction. After having a look at the exponential random number variable, it dawned on me that the argument had to be 1/mu, instead of mu.

The exponential random variable is computed to know when an event actually happens. Apparently, most of the time nothing will happen.The issue with the number of simulations was solved by lowering the initial number of susceptible individuals, and increasing the m parameter.

This is a screenshot of the result:

This is the worksheet: MaplePrimes_SIR_model_simulation_Gillespie_algorithm.mw

I guess the result is correct now: the infected community decreases down to zero.

I also started doing a version with procedures where there is less code in the loop. For this version, it looks like the Statistics package has some issue when trying to compute the addition of the list returned by one of the procedures.

Although I've read that the Statistics package needs to work with floats, the commands do work outside the for loop. I'm also using the hfloat option in the loop, so I don't know what's the issue with this.

This is the attempt: MaplePrimes_SIR_model_simulation_Gillespie_algorithm_proc.mw

Any ideas?

Thanks,

jon

@Carl Love 

Dear Carl,thanks, I did not know that restart ahd that effect on the random number generator.

That definitely makes the program yield different results.

However, I still need to test more thant 2000 time steps to get some sort of meaningful trend in the plot, and it should be enough with a 0..100 range.

Any guess?

This is the updated worksheet: MaplePrimes_SIR_model_simulation_Gillespie_algorithm.mw

Thanks,

jon

@Carl Love 

Carl,

thanks again for your endeavor. I think you were quite right.

However, there must be some other error: the process is not stochastic: I always get the same values, the same plot. I should be getting different realizations each time. Although all would follow the same trend, they should lookdifferent, as in the picture attached a couple of messages ago.

Any ideas?

Besides, I should be getting the trends much earlier than what I actually do: I am forced to set an upper limit of tmax:2500 to see some meaningful plot. That number should be as low as 100. By tmax:=100, and i=30 the infected poupulation (In) should be approximately zero.

This is the current plot (the horizontal axis is limited to the last value of the counter of the loop):

This is the updated worksheet: MaplePrimes_SIR_model_simulation_Gillespie_algorithm.mw

Thanks,

jon

@sazroberson

Just in case there was a misunderstanding, the pictures attached to the previous message were not the result of the Maple program. They are the results I should be getting. But I'm not.

Any ideas why the simulation is not truly stochastic although I'm using random numbers?

Thanks,

jon

@Carl Love 

Hi there,

any updates on this, please?

Thanks,

jon

Hi there,

any updates on this, please?

thanks,

jon

@sazroberson 

Well, I was thinking in ai, which what Kitonum's answer gets, if I'm not mistaken.

But if gettin ic is easy enough, then for the sake of completion, it would be useful to post it for the rest of the mapleprimers.

Thanks,

jon

@Kitonum 

Dear Kitonum,

OK, that worked.

And yes

Search(2,L)

in this case would yield

2

the position of the "2" element in the L list.

I think that will make it.

Thanks,

jon

@sazroberson 

Sorry, sazroberson, I forgot to specify; they are numeric elements.

@sazroberson 

Dear sazroberson,

yes, your are right; I was not propagating the values.

However, it looks like it is just a part of the problem: the simulation does not get the shapes for S, In and Rec I should be getting. Have a look a this picture: the left-hand is for 100 simulations; the right-hand is the mean.

This is the new worksheet: MaplePrimes_SIR_model_simulation_Gillespie_algorithm.mw

Thanks again,

jon