@itsme You're right - I'll try to add it for a future version of Maple.

Erik.

@Mac Dude

Yes, in order for the code to work correctly you have to provide all inflection points. The code should not care about additional points (actually, any additional points, whether or not the 2nd derivative is 0).

I think you're correct that the code is less efficient than a simpler solution if the PDF changes for every invocation. If efficiency is very important, it might be possible to use some of the ideas there to get a tight fitting envelope in your situation, though I think it depends on being able to describe the inflection points relatively easily.

Recalling writing those posts, I think I just didn't consider non-normalized PDFs. The examples I had in mind were all easy to normalize so it just didn't come up. (In the context of the Statistics package, you would definitely need to normalize the PDF - there, everything relies on it being normalized.)

With regards to your PS: I don't know of an easy way to find such posts. I expect most of them would be "featured posts" when they are written, but there's no easy way to find them after the fact. Another favourite of mine is Joe Riel's post about Cartesian products of lists, at http://www.mapleprimes.com/posts/41838-Cartesian-Products-Of-Lists. I reread it once every year or so :)

Erik Postma
Maplesoft.

Hi Mac Dude,

First off, I'm glad to hear you enjoyed my post. Here are a couple of notes in reply to your ideas and suggestions:

1. The process you arrived at looks sound, pretty efficient, and correct to me. That's what matters most!
2. Most of the time spent in generating the 600 points you show, is undoubtedly spent in finding the inflection points of the PDF. This is a step that's needed to make the particular implementation of Maple's envelope generation work. This implementation is explained in a series of four mapleprimes posts I posted four years ago. Finding the inflection points is done in pure Maple code, and it can be slow for complicated PDFs; once that's done, the actual generation of points is done in external C code, it requires very few PDF evaluations, and it's quite fast.
3. The technique from the end of the Sample help page where one leaves a parameter unassigned when calling Sample, and then sets it before calling the resulting procedure, is not efficient for use with the envelope rejection generator: the envelope needs to be recomputed for each invocation anyway, so there's no noticeable benefit to having the sampling procedure pre-computed. (For builtin distributions, it is efficient, because those have pre-written C code for generation of random numbers without substantial set up, and it's just a matter of sending the correct value in.)
4. If you know something about the conditional PDF for y given the value of x, that tells you where one can find its inflection points (as a function of x), then you can adapt the code given in part four of the series linked above to possibly speed the random generation up a bit more. As always, the question is, how much do you need the extra speed vs. how much time would it cost to write this adapted version.

Let me know if you have more questions or remarks.

Erik Postma
Maplesoft.

Thanks for your question, @Majmaj.

The reason that Maple uses the formula that it uses, is the same reason that it divides by N-1, rather than by N, to obtain the variance: this is the formula that gives you an unbiased estimate of the distribution variance -- at least, if the distribution that your values come from is normal. This is really the only formula that makes sense for the general case: where the weights can be any positive numbers. (Consider, for example, if your weights are 1/10, 1/9, 1/8; then using sum(w[i])-1 does not make sense.)

You could argue, as @CarlLove does below, that Maple should detect the case where all weights are integers (or the weights sum to a number greater than 1, or some other condition) and in that case use a different formula - the formula that assumes that weights are just repetitions. I don't agree with this proposal - I think using formulas with different outcomes, in a manner that is controlled by the values of the arguments, and not even in a continuous manner, is just bad practice; it would make the behaviour of Maple much less predictable than it is currently.

If you want to view weights as repetitions in the Statistics package, then it's best to do the expansion yourself: if your data set and the numbers of repetitions are both given as lists, you can do this nicely as follows:

L := [3,2,2];
W := [7,5,8];
LW := zip(`$`, L, W);
Statistics:-StandardDeviation(LW); # returns 0.489360...

Hope this helps,

Erik Postma
Maplesoft.

Note that both ?copy and ?LinearAlgebra:-Copy will work. They do almost exactly the same thing if their argument is a Matrix or Vector.

Erik Postma
Maplesoft.

Note that both ?copy and ?LinearAlgebra:-Copy will work. They do almost exactly the same thing if their argument is a Matrix or Vector.

Erik Postma
Maplesoft.

Hi herclau,

If you call ?randomize with a fixed argument, you always set the random generator to the same state, so you will always get the same matrices. That's its purpose.

If you want to get different results, it's probably best to just not restart in between subsequent runs. Alternatively, you can omit the argument in the call to randomize(); that will set the random seed to something based on the current internal clock time of your computer (though I believe it only uses a granularity of seconds, so don't do it in, say, a loop - in which case you're doing something wrong anyway, because for obvious reasons you can't restart from within a loop.)

HTH,

Erik Postma
Maplesoft.

Hi herclau,

If you call ?randomize with a fixed argument, you always set the random generator to the same state, so you will always get the same matrices. That's its purpose.

If you want to get different results, it's probably best to just not restart in between subsequent runs. Alternatively, you can omit the argument in the call to randomize(); that will set the random seed to something based on the current internal clock time of your computer (though I believe it only uses a granularity of seconds, so don't do it in, say, a loop - in which case you're doing something wrong anyway, because for obvious reasons you can't restart from within a loop.)

HTH,

Erik Postma
Maplesoft.

Hi Daniel / SamuelTuvare,

Could you show us the details of what you have done so far?

Erik Postma
Maplesoft.

Taking the time interval between t=0 and t=0 and dividing it into 100 time steps, you get steps of length 0. The system is not set up for this. Similarly, you also can't ask for the value at negative times.

Hope this helps,

Erik Postma
Maplesoft.

Taking the time interval between t=0 and t=0 and dividing it into 100 time steps, you get steps of length 0. The system is not set up for this. Similarly, you also can't ask for the value at negative times.

Hope this helps,

Erik Postma
Maplesoft.

@Markiyan Hirnyk : Another option that is equally fast (measured to within measurement accuracy on an i7 running linux-64) is:

nops(select(`<`, thelist, 0));

The big advantage of this over the version with c -> c < 0 is that it only uses a kernel-internal procedure in the inner loop. 

Met vriendelijke groet,

Erik Postma
Maplesoft. 

@Markiyan Hirnyk : Another option that is equally fast (measured to within measurement accuracy on an i7 running linux-64) is:

nops(select(`<`, thelist, 0));

The big advantage of this over the version with c -> c < 0 is that it only uses a kernel-internal procedure in the inner loop. 

Met vriendelijke groet,

Erik Postma
Maplesoft. 

 It uses default values of 10^2 (equally long) time steps and 10^4 iterations.

This is determined by the procedures Finance:-ExpectedValue:-ProcessParameters and Finance:-ExpectedValue:-ProcessCommonOptions. You can look at the source code for these guys by doing:

kernelopts(opaquemodules=false):
showstat(Finance:-ExpectedValue:-ProcessParameters);
showstat(Finance:-ExpectedValue:-ProcessCommonOptions);

and if you showstat(Finance:-ExpectedValue) you can see how these two are called.

Hope this helps,

Erik Postma
Maplesoft.

 It uses default values of 10^2 (equally long) time steps and 10^4 iterations.

This is determined by the procedures Finance:-ExpectedValue:-ProcessParameters and Finance:-ExpectedValue:-ProcessCommonOptions. You can look at the source code for these guys by doing:

kernelopts(opaquemodules=false):
showstat(Finance:-ExpectedValue:-ProcessParameters);
showstat(Finance:-ExpectedValue:-ProcessCommonOptions);

and if you showstat(Finance:-ExpectedValue) you can see how these two are called.

Hope this helps,

Erik Postma
Maplesoft.