@mmcdara this is better than any textbook I own and have read on probability distribution with Maple. I can't thank you enough.

Yes, unfortunately your statement:  

"But the problem is that Maple meets difficulties to compute some ot these statistics if not returning wrong results (I'm almost sure of that but I need to check it one more time)."

Is true. I tried with a several tests, the theoretical Median and most Quantiles returned complex numbers, while the outcome of the sampled (empirical) 10^5 are consistently real numbers.

Anyways. Appreciate very much the help @mmcdara . I have homework to catch up, All_demonstrations.mw being one of them. "See you" much later.

Thank you very much @mmcdara . This helps and solved it. New things I learnt from this exercise that I would not get (or very hard to find) are:

1. It is syntax-ly (sorry for the odd term) correct in the use of 2 stages of unapply function [in the local variable declaration and inside the Distribution function - ie. local pdf := unapply(PDF......) and PDF = unapply(pdf....)]. Cool

2. To apply the location (shift) and scale parameters. The key is to multiply (transform) the underlying distribution with (t-a)/b . Where a=location, b=scale.  I thought I saw this in my first question in this forum, but did not quickly get it.

3. To define the Variance inside the Distribution function, we have to square the scale (like b^2) and "multiply" with the underlying distribution.

4. There are no transformation factors applied to Skewness and Kurtosis inside the Distribution function

5. The Mean and Median inside the Distribution function are the least difficult to understand. ie. just add the location to "shift" the results after applying the scale.

Now, I was wondering if there is a single book reference (maybe a mathematical statistics book? Or a distribution function transformations book?) that explains more on 2., 3., and 4., and 5. above? 

Thanks again, and double awesome.

Thank you very much. I will experiment more on this and will sure ask more questions in due time on another thread/post. With two question threads in this week alone, I never thought I get this much and thorough help. Awesome, @mmcdara

@mmcdara I did not know that default values can be defined in the arguments of a procedure, thanks for this tip.

Being new to Maple's procedures, I was thinking if I could debug a procedure easily. If one wants to know the output of every local declarations like removing the colon (:) after the end of each command line inside a procedure?

Let me know if this is a new thread I could post? I am still trying to understand the TriangularTypeII procedure (before implementing to an exponential distribution (with the same 2 pairs of quantiles and probabilities)). 

@mmcdara I think I got it. To understand fully, let me try with local eqs := {F(Q[1])=Q[2], F(P[1])=P[2]} , to confirm my understanding.  I very much appreciate your help @mmcdara

This makes the solving of a and b foolproof @mmcdara , plus I get to learn the map function from a real case.  

Sorry for my slow logic, but I still dont get (in the original TriangularDistribution.mw and the version 2) how Maple realize that the pairs of Q and P are quantiles and their probabilities (the 2 pairs) in this equation:  

local eqs := {F(Q[1])=P[1], F(Q[2])=P[2]}

Because Q and P could be anything as far as the software is concerned?  The bunch which leads to the above equation are:

   local TRI := RandomVariable(Triangular(a, b, c)):
   local F   := unapply(CDF(TRI, t), t):

How does Maple know that Q and P are pairs of quantile and its probability? In other words, how does this get defined so Maple knows this?  

Beautiful solution, @mmcdara . Appreciate it.  As basic as it may seem, I am just made aware of the powers of the "op" function. ... And I pondered on these statements, which I believe is key to the reparameterization:  

   local F   := unapply(CDF(Triangular(a, b, c), t), t):
   local eqs := {F(Q[1])=P[1], F(Q[2])=P[2]}:
   local AB  := solve(eqs, {a, b}):
   local T   := RandomVariable(Triangular(op(eval([a, b], AB)), c)):

Yes, on the challenges to get "a" and "b" numerically with 5 knowns, I got 4 pairs of solutions with "solve",

{a = 0.5857864376, b = 7.414213562}, {a = 3.414213562, b = 4.585786438}, {a = 1.784749563, b = 4.451416230}, {a = 3.548583770, b = 6.215250437}

but I was puzzled that I got the right pair of solution using "fsolve".

One last question on this thread, is there a Maple documentation on how this got defined? 

local eqs := {F(Q[1])=P[1], F(Q[2])=P[2]}:

Thank you much.

This has been a very thorough approach on the subject of Beta reparameterization, and I can't thank you enough, @mmcdara .

I particularly learnt (was initialy lost, but found out by searching some subject in mathematical statistics) was this steps:

MO := Mode(U) assuming alpha > 1, beta > 1;

reparam := solve({ME=MyMean, MO=MyMode}, {alpha, beta}); 
eval(reparam, {MyMean=(me-P)/(Q-P), MyMode=(mo-P)/(Q-P)}):
simplify(%);

MyMean=(me-P)/(Q-P) and MyMode=(mo-P/(Q-P)

which are key to reparameterization and wonder if I could ever found this on my own (up to the correct Maple syntax). 

The addition to End-point estimators is also a bonus that I am sure that I will bump into in the near future.  It is relevant to what I am going to do next, which is recreating these (just like the Beta reparameterization): Vose and @Risk , which evalauates numerically the minimum and maximum based on the given low and high percentiles.  Wish me luck on this steep learning curve, thank you very much once more @mmcdara .

Appreciate the patience and perseverance in tackling this @mmcdara. As my first post in the forum I did not know how to respond to the threads, and have no idea how many will participate to contribute. BSD_2b.mw did solve this and taught me several lessons that is very difficult to find in the web in the Maple language.

@mmcdara by the way, the 3rd and 4th central moments in file Use_formulas.mw (skewness and kurtosis) should have been 0.14606 and 2.02666 respectively and checked by the Vose's and Palisade's website of the skewness and kurtosis formulas. Here succesfully recreated in Maple: 

Download Distribution-Beta-Subjective3.mw

@mmcdara sorry to drag on with the subject, but the ultimate goal of being able to simulate from a custom defined distribution of SB(minimum, mode, mean, maximum) has not yet been achieved.

Let's point out to the file you proposed in the previous thread: Use_formulas.mw. The file stopped at calculating the 25 percentile of the created analytical distribution, which is plotted in plot(f, x=Min..Max, gridlines=true). A successful Distribution function (Distribution(T1, T2, .., Tn)) should enable the use of the function Sample(X, n, numeric_option, output_option) . I have referred to the Maple help file and took the excerpt of it below. Please take note of the yellow highlighted comment in the mw file. By the end of the day (ideally) we should be able to plot a Density Plot and an overlay of the sampled histogram of the reparameterized Beta (min, mode, mean, max). Why do I need to go as far as being able to draw random samples of Beta(min, mode, mean, max)? It is to perform Monte Carlo random number generation of this particular Beta(min, mode, mean, max) realizing the possible failures of generating it because of the risks associated with defining this sort of parameters to the Beta as pointed in the threads before this.

Download Custom_Distribution_Template.mw

@mmcdara , I am glad that you found this subject so that I learnt from the methods of solving this that could'be taken weeks for me to approach, let alone solve.

Please dont hate me for what I am going to ask next: on the case of solving the parameters by using the definitions of Mean, Variance, Skewness, Kurtosis, Mode and Quantile, which is the "theoretical" ("by definition") properties, now how do I get random variates of the defined subjective Beta so that the mean is the average of the sampled data etc., approaching the real "theoretical" mean (convergence) when the sample is very large (like _Sample:= Sample(BetaSubj, 10^5); _Average:=Mean(_Sample)? 

Don't think it too much, it may be too much to ask from me. Lol.

Thank you for digging deeper into this, @mmcdara .

Yes, ModelRisk from Vose applied the reparameterized Beta distribution to generate random variates in an Excel add-in for risk analysis. Its competition Palisade @Risk also exactly applies the same definition here: RiskBetaSubj , with a minor typo in the numerator of the alpha definition. So, it's not uncommon that reparameterization (given the risks of failing to calculate) is given to distributions (like the Beta-PERT, for example) so that subject matter experts can easily define their estimates for the purpose of generating random numbers out of that parameters. Most people are more familiar with common concepts of "minimum", "most likely", or "average", and maximum rather than if they are asked to describe an uncertainty using greeks.

Correct, @mmcdara , as per exact pdf definition (and parameter restrictions) of the first thread of this subject (and image copied previously) - with no closed form CDF - so the results are solved numerically (as a monte carlo procedure) with values (min, mode, mean, max) initialized just after restart to enable to evaluate the parameter restrictions.

@mmcdara , appreciate much for your solution. This certainly broaden my perspective and pave a way to the solution I was trying to find. But, just concentrating on WAY 2, the distribution defined is a generalized beta, where the parameters I was really looking for are Min, Mode, Mean, Max (not alpha, beta, min, max), like the link provided in earlier to get exactly these (except for the non closed form CDF). The properties of this subjective beta (not generalized beta) have been applied in 2 Monte Carlo softwares (Palisade's @Risk and Vose - the link above), and I attempted to solve this numericaly (after symbolicaly defining) to replicate what the Monte Carlo softwares were outputing. The mw document I attached checked fine with Min := 3; Mlikely := 8; Avg := 9; Max := 18, but failed with say: Min := 1000; Mlikely := 1400; Avg := 1500; Max := 2100. (numerically as defined in the image, where Mean=mean, Mode=Mlikely, and so forth). So, the B and BSD below need to be adapted for these parameters.

B   := RandomVariable('Beta' (alpha, beta)):
BSD := x__min + (x__max-x__min)*B:

to get parameters of Minimum, Most Likely (Mode), Average (Mean), and Maximum. Rather than alpha,beta, minimum and maximum. Thanks again, @mmcdara !