Test Random Functions

[JamesPHoughton]

It would be good to do an integration test of all of the random functions we support:

• "lognormal": "np.random.lognormal",
• "random normal": "functions.bounded_normal",
• "poisson": "np.random.poisson",
• "exprnd": "np.random.exponential",
• "random uniform": "functions.random_uniform",

The challenge is that our testing infrastructure matches values at each timestep, and that won't work with random functions.

Lack of this sort of testing leads to issues like #105.

[enekomartinmartinez]

Some tests have been added for the supported random functions. We may need to add tests for the different moments and ranges of the distributions to ensure they are being generated correctly.

[hyunjimoon]

Could you add exponential function? My colleagues use it often for lead time and informed me to plead for adding this function :)

Following normal's code from here, and vensim parameterization doc

RANDOM EXPONENTIAL(m,x,h,r,s) provides an exponential distribution starting at 0 with a mean of 1 before being stretched, shifted and truncated.

I tried

"random_exponential": (
"np.random.exponential(scale=%(0)s, size=%(size)s)",
("numpy",)),

[enekomartinmartinez]

Hi @hyunjimoon

We were using np.random.exponential before, but we removed it as it does not do the same thing as Vensim. Vensim stretches, shifts, and truncates most random functions. That's why we use scipy.stats.truncnorm for random_normal instead of other random normal functions from numpy. Maybe the same could be done for random_exponential using scipy.stats.truncexpon. However this class rvs method doesn't allow setting the minimum (m in Vensim documentation a in scipy documentation), then the final distribution may be different is the minimum is greater than the shift.

I could add it, but it will be used at your own risk if you don't mind. We are missing random function tests (that was the purpose of this issue). Ideally, we should have a bunch of data generated by different configurations of the random functions in Vensim. Then, compare with our implementation different order moments (like mean, variance, skewness, and kurtosis) and minimum, maximum, and some percentiles to check that the distribution behaves the same. I currently have no time nor access to Vensim to develop that, so I cannot ensure that it will behave the same.

If you want to test this approach, this is what it should be added in pysd.builders.python.python_functions, note that the first value is ignored as it has no implementation for this function:

    "random_exponential": (
        "stats.truncexpon.rvs(%(1)s, loc=%(2)s, scale=%(3)s, size=%(size)s)",
        ("scipy", "stats")),

[hyunjimoon]

Sounds great - although my collaborator (he told me to pay gratitude to you and @JamesPHoughton for this amazing library) will be the one who'll test, I'll make sure to share whether it works, and if not the error log in this thread.

[enekomartinmartinez]

I guess that as the exponential distribution truncation on the minimum doesn't change the distribution maybe it could be done by setting the loc=np.maximum(%(0)s ,%(2)s)? I am not sure about that but if you can generate distributions and test with different inputs it may it work to reproduce exactly Vensim's behaviour.

This would be possible due to the memorylessness property of exponential random variable.

[parksehoon1971]

Hello everyone?

It's amazing that such great work is being done here. I'm the one who asked to add that function. Thank you for responding quickly. Since I am a modeler, I do not have the ability to modify Python or Vensim functions.

But RANDOM EXPONENTIAL(m,x,h,r,s) in the vensim app provides an exponential distribution starting at 0 with a mean of 1 before being stretched, shifted and truncated. So loc=h and scale=r would be used. The problem is the value of 'b', and this value should be set to 5.3, which is 99.5% of the exponential distribution with a mean of 1. Nevertheless, it is difficult to say that the two functions are exactly the same.

The following picture is the result of comparing the random numbers generated by vensim by parameter. I hope this helps you convert the function.

And I hope this function will be included in the next released version.
thank you!

[enekomartinmartinez]

Hi @parksehoon1971,

I could try to add try to add the function. However, I will be doing it in my free time, so I will need a little collaboration from you. I can develop some tests to compare random distributions., so we can also test the existing ones and potentially add others. I will also add a small toolkit to plot the results so we can visually compare them if something is not working as expected.

The deal is that I will implement the testing tool, and I will ask you for data. We may need to generate a very big bunch of data to test it, and I do not have access to Vensim anymore. If you agree, I will send you some configurations to run (I need to check before which number of points is reasonable so we have strong tests, but they don't take too long). I will ask you to generate data in the same way you did in the example, a TAB file with the output per column from Vensim may work, but will be nice to have in the first row the function that was called (e.g. RANDOM EXPONENTIAL(0, 100, 0, 1, x) ).

[enekomartinmartinez]

I started working on the test at https://github.com/SDXorg/pysd/tree/random_functions
Given scipy limitations and GitHub file size limitations I guess 10e5 is a good number of samples, I have been testing it for the same sample mean, variance, skewness and kurtosis. For all of them, the tests pass with 1e-5 absolute and relative error.

[parksehoon1971]

Hi @parksehoon1971,

I could try to add try to add the function. However, I will be doing it in my free time, so I will need a little collaboration from you. I can develop some tests to compare random distributions., so we can also test the existing ones and potentially add others. I will also add a small toolkit to plot the results so we can visually compare them if something is not working as expected.

The deal is that I will implement the testing tool, and I will ask you for data. We may need to generate a very big bunch of data to test it, and I do not have access to Vensim anymore. If you agree, I will send you some configurations to run (I need to check before which number of points is reasonable so we have strong tests, but they don't take too long). I will ask you to generate data in the same way you did in the example, a TAB file with the output per column from Vensim may work, but will be nice to have in the first row the function that was called (e.g. RANDOM EXPONENTIAL(0, 100, 0, 1, x) ).

@enekomartinmartinez

Thank you your efforts!
And I'm sorry to hear that you use up your free time to work. If I could do what you are asking, I would.
I would like to know specifically what task you are requesting. Is it okay for me to create a tab file with the results for the parameters you set? Or can I just decide the parameters arbitrarily and create random numbers?

[enekomartinmartinez]

Hi, @parksehoon1971

I have been doing some other tests. I guess the better approach is to compare both with Vensim but using a bigger tolerance and then with expected values calculated analytically so we can test with less error.

In this process, I have discovered that the current implementation of RANDOM NORMAL may not be working as in Vensim, but I need to test against Vensim data first; I think it should be translated as stats. transform.rvs((%(0)s-%(2)s)/%(3)s, (%(1)s-%(2)s)/%(3)s, loc=%(2)s, scale=%(3)s, size=%(size)s) instead of stats.truncnorm.rvs(%(0)s, %(1)s, loc=%(2)s, scale=%(3)s, size=%(size)s). This is because scipy does the truncation at the beginning and Vensim at the end. I guess the same will happen with truncexp.

I will ask you to produce the following files with the first column Time and 4 more columns (1 per function) total 5e5 data (you can just use 500000 as a FINAL TIME 0 as START TIME and TIME STEP = SAVEPER = 1:

uniform_vensim.tab

RANDOM 0 1()
RANDOM UNIFORM(0, 1, 0)
RANDOM UNIFORM(-3, 10, 0)
RANDOM UNIFORM(0.5, 7.8, 0)

normal_vensim.tab

RANDOM NORMAL(-10000, 10000, 0, 1, 0)
RANDOM NORMAL(0, 1000, 0, 1, 0)
RANDOM NORMAL(-100, 100, -3, 0.5, 0)
RANDOM NORMAL(3, 6.9, 5, 5, 0)

exponential_vensim.tab

RANDOM EXPONENTIAL(0, 10000, 0, 1, 0)
RANDOM EXPONENTIAL(2, 100, 0, 1, 0)
RANDOM EXPONENTIAL(0, 100, 2, 0.5, 0)
RANDOM EXPONENTIAL(0, 7.2, 5, 5, 0)

If you have doubts maybe is easier to talk via slack, you can join the workspace we have sdtoolsandmet-slj3251.slack.com

Files should be something like:

TIME	RANDOM NORMAL(-10000, 10000, 0, 1, 0)	RANDOM NORMAL(0, 1000, 0, 1, 0)	RANDOM NORMAL(-100, 100, -3, 0.5, 0)	RANDOM NORMAL(3, 6.9, 5, 5, 0)
0	-0.5410708424521473	1.2343302669140448	-2.908936457110545	5.371382640882544
1	-0.10606290064934429	1.3294622453972802	-3.410060359110077	4.590287428710126
2	-0.22508427875464468	1.081605702435474	-3.2322681729925558	6.71351847321953
3	0.4878344429732129	0.9297296474845612	-3.0218147585434414	4.350613094168536
4	0.05285961565049486	0.017638019804912756	-2.7294256353078876	5.4415228651410334
5	-0.1622428174579056	1.9809308250249869	-3.190120157158069	4.965529602645091
6	0.5523188674122674	0.47347386147119175	-2.511506059782375	3.8237159295221783
7	0.06799519255688485	0.12374653822080052	-2.9745080657171785	5.776546033861598
8	-0.8434038278440443	0.5133261881921402	-2.805638055759695	4.7784330939969
9	1.0499652568414048	0.5808988226499702	-2.8817495450573634	3.3430309367958433
10	0.4109332670053265	1.070439555612693	-3.1794599928504805	4.566733281204497
11	-0.35158163268601045	0.22292757847150244	-2.290255558772565	5.577255856832461
12	-1.2769363011471542	0.5482801198168058	-2.800530300045481	5.7131187035792355
13	0.8295810596274037	0.14530007466461276	-3.158529102136417	3.120575688911325
14	1.2439083749334197	0.14307824759338328	-3.4669227457646605	5.35811156139726
15	-0.1792073396031722	0.05228387900947716	-2.235824330658638	5.609613372287555
16	0.05868418926452117	0.04444056128945642	-1.695134914263651	5.925529807803697
17	0.3814385609858069	0.9147874534708827	-2.276413092448389	5.476336168668346
18	0.9233599484369635	0.9375131240626801	-2.589213082613727	3.1081070731028118
19	-2.084790950770809	1.0973204395399871	-2.933329448666583	4.0241614880449665
20	0.5235736299894298	1.3052337668345655	-2.7861911250840876	5.532444816898912
21	1.7396153692964915	0.02683788205839089	-2.782761700036949	6.25835846278721
22	-1.1412124568777826	1.3237701767239298	-2.676355811554529	6.440757641439723
23	0.8352594376097893	0.3735815809880256	-3.1780540890024946	5.3872484893237464
24	0.013082723511630121	0.7868972334868717	-3.3406568442115447	6.305734490503385
[...]

[parksehoon1971]

Hi, @enekomartinmartinez

I'm sorry for the delay in communication due to the time difference. I have finished what you asked me to do. Please test using the attached file. thanks!

tab_files.zip

[enekomartinmartinez]

Hi, @parksehoon1971, as soon as the tests pass, I will create the new release. Hope it helps.

Just remember that if you publish anything using PySD you can always cite our paper

[parksehoon1971]

@enekomartinmartinez, Thank you so much for updating the library so quickly. This update removed one obstacle in implementing the web app we were trying to develop. I hope you are always happy.

[parksehoon1971]

@enekomartinmartinez,
Thanks to you, I am able to implement the Vensim model on the web. However, when random numbers occur, the seed value keeps changing, so the same result cannot be produced again. It would be of great help if the seed value of the random number could be fixed. Please refer to the attached two pictures showing different simulation results in the same settings. I'm sorry to ask you a favor when I know you're busy. Still, I hope you have a nice weekend. thank you!

[enekomartinmartinez]

Hi @parksehoon1971

try using numpy.random.seed before running the model, this should set the seed for all the calls to random functions.

[parksehoon1971]

@enekomartinmartinez

Yes, thank you for helping me solve the random number problem.

[parksehoon1971]

Hi, @enekomartinmartinez

While running the model with PySD, something strange occurred. There was a difference between the data executed by Random Normal() in Vensim and the data executed in PySD. In particular, there is a large difference between the upper and lower limits.

mean daily demand = 100
customer order = RANDOM NORMAL(0, mean daily demand2, mean daily demand, mean daily demand2/3, seed)

So I created only one variable and generated a random number, and there was no problem.

I've attached the file with the problem, so if you have the time, I hope you can solve this problem. Thank you for always!

q_system.zip