Although the worksheets provided here have been developed under Maple 2015, they should work correctly with newer versions, except perhaps for commands that use the 'op' function ('piecewise' mainly).
In the sequel the acronym 'pdf' stands for 'probability density function'.
CONTEXT
This post originates from a recent question by @JoyDivisionMan and the ensuing discussion.
In a few words, the OP noticed Maple 2025 failed to return a result and asked why. In his reply, @acer identified a code regression somewhere in between Maple 2023 and Maple 2025.
Indeed Maple 2023 (which I do not have) provided an answer, just so "my" Maple 2015 did... but thee problem was that Maple 2015's answer was wrong (one might aske what is worse: pitifully admit one can not find an answer or provide a wrong answer?).
It is likely that Maple versions prior to 2024 both return an incorrect answer given that questions relative to Statistics, and PDF in particular, are so marginal on this site that it is highly improbable that a Software Change Request has ever been issued on the subject.
In this post I explain how we can calculate the result by hand (only elementary maths required), why Maple 2015 (and likely newer versions) returns an incorrect result, why Maple generally fails in returning a result, and finally provide several examples to illustrate that even mathematically simple, the construction of the pdf often requires human intervention (unless a sophisticated AI-driven algorithm is used... at lest in my opinion).
To be clearer, let X some continuous univariate random variable (CURV) and 𝟇 a real valued function from ℝ to ℝ (in the most general sense). Let Y the random variable defined by Y = 𝟇(X).
Here is the main claim: unless very specific situations, you must consider that Maple cannot correctly build the pdf of 𝟇(X) if 𝟇 is not a monotone function (and even sometimes a strictly monotone function).
A few noticeable exceptions are
X any CURV , 𝟇 : x ⟼ x2 (correct solution even if 𝟇 is not monotone)
X ~ Uniform(-1, 1) , 𝟇 : x ⟼ arctanh(x) (no result returned even if 𝟇 is strictly monotone).
In @JoyDivisionMan's original question it is only by chance that Maple provided the correct result for the non monotone function 𝟇 = cosine ... but it failed for sine function, which motivated the question.
A GOOD DRAWING WORTH A THOUSAND WORDS
Here is a picture to help understand how to get the pdf of Y = 𝟇(X) for a non monotone function 𝟇 (the case of a monotone function directly comes from this latter).
In this illustration X ~ Uniform(0, 2𝜋) , 𝟇 : x ⟼ sine(x).
To ease the explanation, I write X as a mixture of three uniform random variables X1, X2, X3, whose supports are the intervals of the three branches of 𝟇. More formally, X = (1/4)∙X1 + (1/2)∙X2 + (1/4)∙X3.
The restrictions of 𝟇 to these three branches are denoted 𝟇1, 𝟇2, 𝟇3.
The large rectangles below the horizontal axis represent the pdf of X1, X2, X3 and the blue curve the 𝟇 function.
The image of the interval [y-dy, y+dy] by the inverse functions 𝟇1(-1) and 𝟇2(-1) of 𝟇1 and 𝟇2 are represented by the vertical rectangles [x1-dy, x1+dy] and [x2-dy, x2+dy] .
These two intervals bring a contribution to the pdf of Y (light gray blue on the right) represented by the horizontal violet rectangle on the right side of the picture.
The probability Prob(Y ∊ [y-dy, y+dy]) that Y belongs to the interval [y-dy, y+dy]) is simply the sum
Prob(Y ∊ [y-dy, y+dy]) = Prob(X1 ∊ 𝟇1(-1)([y-dy, y+dy])) + Prob(X2 ∊ 𝟇2(-1)([y-dy, y+dy]))
Let 𝟇'b(-1)(y) denote the derivative of 𝟇b(-1)(y).
Making dy tends to 0 gives
pdf(Y = y) = pdf(X1 = 𝟇1(-1)(y)) × | 𝟇'1(-1)(y) | + pdf(X2 = 𝟇2(-1)(y)) × | 𝟇'2(-1)(y) |
As I said before there is truly no big math behind this, Except maybe those absolute values?
To understand where they come from zoom in on the rectangle 𝜔 = [x1-dx, x1+dx] ╳ [y-dy, y+dy] and denote X𝜔 and Y𝜔 the restrictions of X1 and Y to 𝜔.
Locally Y𝜔 is proportional to A+B∙X𝜔 where constant B = 𝟇'(x1) and the value of constant A does not matter here.
So the pdf of Y𝜔 is (a classical result) pdf(Y𝜔 = y) = pdf(X𝜔 = (y-A)/C) / |C|.
Few details can be found Here.
MAPLE FAILURES AND WEAKNESSES
So why did Maple, at last some versions, produce a wrong result and why some versions are not even capable to return one?
The reason is that there is no big math only at first sight...because determining the inverse function of 𝟇 can be quite tricky as soon as 𝟇 is not one-to-one map, for instance when 𝟇 is not strictly monotone.
When it is so 𝟇(-1) must be defined for all the branches whose definition intervals intersect the support of X.
I spent a lot of time debugging the procedure Statistics:-PDF to understand why it either fails or produces incorrec results.
The sine_debug_nodebugoutput.mw worksheet presents the "X ~ Uniform(0, 2𝜋) , 𝟇 : x ⟼ sine(x)" case (as Mapleprimes stubbornly refuses to upload the worksheet containing the debugger trace, I convert it to sine_debug.pdf to help you see this trace).
To orient the core development team correcting this procedure (assumming they care), the critical procedures are
Statistics:-RandomVariables:-PDF:-Univariate:-GetValueTab[anything]
and Statistics:-RandomVariables:-GetInverse
At the very end it is this second procedure which is truly responsible of Maple failing to provide a result or returning an incorrect on, because it does not correctly build the inverse functions 𝟇b(-1) for all the branches b which matter.
I wrote above that "it is only by chance that Maple provided the correct result for the non monotone function 𝟇 = cosine". Indeed Statistics:-PDF returns a wrong result when X ~ Uniform(z, z+2𝜋) and z is not a multiple of 𝜋/2 (see cosine.mw).
Other important situations where Maple fails returning a result are those where 𝟇 is a polynomial function with different zeros located in the support of X.
I did not trace them but it seems that Statistics:-PDF does not know how to build the 𝟇b(-1) in this case (even though it is quite simple, see "Polynomial" examples below).
A SELECTION OF EXAMPLES
Here is a selection of examples to demonstrate that even in rather complex cases the pdfs of 𝟇(X) can be constructed quite easily (note that Maple either fails to compute them or to provide a correct result):
OPEN QUESTION
Tracing Statistics:-PDF reveals an already complex algorithm designed to handle a broad variety of "canonical" situations. Even this the algorithm fails in almost all non-toy-problem such as this compilation proves Maple_failures.mw.
At first sight, there seems to be a contradiction between the (apparent?) simplicity with which one can obtain, by hand in sometimes in an ad hoc way, the expression of pdf(𝟇(X)), whether it is exact or truncated, and the complexity of the Statistics:-PDF algorithm, which results in failure in all non trivial cases.
This observation leads to the important question "Is it possible to rewrite Statistics:-PDF in order to enlarge its domain of success?".
I have the feeling that this means designing an algorithm which focuses more on mimicing the human reasoning than identifying "canonical" situations (as it is done today).
An AI-driven algorithm maybe?
I repeat here that the maths are very simple, and all the more simple if you represent the random variable X as a mixture of components X1, ... XB, both having the same (truncated) distribution than X, and whose supports identify to the intervals of definition of the B branches of 𝟇 over the whole support of X.
The only difficulty lies in the identification of these branches and in the construction of the functions 𝟇b(-1) over the supports of each Xb.
I have no answer to this question.
