@NanaP 

I'm not sure I understood correctly so correct me if I'm wrong:

• You have a random variable X with a SB (short for "subjective Beta") distribution of parameters xmin, xmax, alpha, beta (I'm not sure of this point, maybe they are xmin, xmax, mode, mean ?).
• You have a sample S of X, let's say of size N, from which you compute some empirical statistics, for instance
  the mean 

  m := (1/N)*add(s[n], n=1..N)
  

• Some facts are:
  • the above empirical statistics is the best unbiased asymptotic estimator of the Mean(X) (of course you already knew that);
  • the classical estimation of the variance 

    v := (1/(N-1))*add((s[n]-m)^2, n=1..N)

    also has this same property.

  • The case of Skewness and Kurtosis is a little bit more complex. If I'm not mistaken there are no formulas for unbiased estimator (nearly unbiased estimators do exist).

  • The case of xmin and xmax is even more complex.
    I have some material concerning what is called "the endpoint distribution", that is the estimation of xmin or xmax from sample S.
    Here is a first reference, https://arxiv.org/pdf/1412.3972.pdf, do not hesitate to ask for some more.

  • Concerning the estimation of the mode there is no unbiased estimator, but a lot of estimators do exist.

@NanaP 

Thank you for introducing me to something I didn't know about, even though I've been working in a parallel domain for a long time.

Hope my contribution will be useful

@acer 
I had observed on my own that isolate didn't work for this example.
So you are right saying that the imaginary unit is not the cause of my difficulty in the sense that I could faced them without having I.

It's just that for the particular case b=1 isolate works if c <> I, and not instead, so my quick jump to an erroneous conclusion.

@acer 

Thank you acer for all these informations.

@NanaP 

The case 

Min := 1000; Mlikely := 1400; Avg := 1500; Max := 2100;

can be solved by using the definitions of Mean, Variance, Skewness, Kurtosis, Mode and Quantile:
 

Download Use_formulas.mw

Note that the mode is not necesarily  unique (Beta(a, a) has two modes at x=0 and x=1 when a < 1 and a continuum of modes whan a=1).
So, I don't know what you want to do, but this re-parameterization you want to use is extremely dangerous.

Download Distribution-Beta-Subjective_Failure.mw

About the Beta distribution there are only, up to my knowledge, these parameterizations:

• shape (a) and scale (b), the most common ant the one Maple adopts;
• mean (M) and variance (V);
• mean and precision (inverse of variance, noted tau below);

As you can see there is a 1-to-1 map between each of pair of parameterization

m := Mean(Beta(a, b)):
v := Variance(Beta(a, b)):
solve({m=M, v=V}, {a, b});

solve({m=M, 1/v=tau}, {a, b});
        /        / 2        \      / 2        \        \ 
        |      M \M  - M + V/      \M  - M + V/ (M - 1)| 
       < a = - --------------, b = -------------------- >
        |            V                      V          | 
        \                                              / 
  /      3        2              / 2                \        \ 
 { a = -M  tau + M  tau - M, b = \M  tau - M tau + 1/ (M - 1) }
  \                                                          / 

This VOSE article you refer about seems to be oriented to risk assesment, for which this Subjective Beta Distribition might be, according to the authors, and likely under some restrictions an interesting tool to achieve this goal.
I'm working on the domain of Reliability Analysis (assessment of probability of failures due to extreme events in engineering, something probably not that far from the "risk assessment" VOSE talk about?]) for more than 25 years and I had never heard of the approach developped by VOSE (but each domain develops it's own method with little inter-operability).

Some intersting links:

Probability ontology
http://www.math.wm.edu/~leemis/2008amstat.pdf
https://en.wikipedia.org/wiki/Relationships_among_probability_distributions

Different distributions
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5013898/bin/supp_btw170_Swat_et_al_Supplementary1_ProbOnto_KBprintout.pdf

Different parameterizations
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5013898/bin/supp_btw170_Swat_et_al_Supplementary2_ProbOnto_useInPharmML.pdf

@NanaP 

Does it mean you want to re-parameterize a Generalized Beta distribution by min, max, mode, min instead of alpha, beta, min and max ?

If it so you will have a problem because the relation (alpha, beta) --> (mode, mean) is not a 1-to-1 mapping.
Take the simplest case alpha > 1, beta=alpha, xmin=0, xmax=1 (thus the Generalized Beta is just Beta(alpha, alpha)).
As Mean(Beta(alpha, alpha)) = Mode(Beta(alpha, alpha)) = 1/2 whatever alpha > 1, a re-parametrization of the Beta distribution in term of mode and mean doesn't say anything about the pdf, cdf, variance, ... of this distribution.

@C_R 

but the expression 

I*Int(f(x), x) = something - 2*I*Int(f(x), x);

comes from some computation and I didn't want to extract the Int component from it's left hand side and isolate it from the equality.
In fact I wrote f(x) but I don't even know a priori what f(x) is.
So you can do something like this:

expr := I*Int(f(x), x) = something - 2*I*Int(f(x), x):
J := select(has, [op(lhs(expr))], Int)[]:
isolate(expr, J)
                    /                         
                   |             1            
                   |  f(x) dx = -- I something
                   |             3            
                  /                           

The problem is that you replace I by some constant c, then isolate works as expected:
 

c*Int(f(x), x) = something - 2*c*Int(f(x),x):
isolate(%, lhs(%))
                   /  /        \              
                   | |         |   1          
                 c | |  f(x) dx| = - something
                   | |         |   3          
                   \/          /              
I*Int(f(x), x) = something - 2*I*Int(f(x),x):
isolate(%, lhs(%))
          /  /        \                   /  /        \
          | |         |                   | |         |
        I | |  f(x) dx| = something - 2 I | |  f(x) dx|
          | |         |                   | |         |
          \/          /                   \/          /

@sursumCorda 

Thank you for your explanation.

https://www.mapleprimes.com/questions/236509-How-Do-I-Solve-This-One.

As I'm in a good mood today, here's a solution

Download Maple_expression.mw

Can't you provide a Maple code for us to see what you are capable of and then guide you towards the solution?
The formulation of your question sounds like "do my homework, I don't want to bother with it".

@KIRAN SAJJAN 

Could you explain more precisely what you want because your last worksheet doesn't contain any error and your request is not clear (at least for me).
I imagine M is the Mach number, Cf is the friction coefficient and Nu is the Nusselt number, aren't they?
Do you mean that you want to plot , for instance, aNusselt number defined by 

(-a4-Rd)*diff(T(Y), Y)

for different values of the Mach number?

Could it be this:
 

Download AGM-2_mmcdara_1.mw

@ogunmiloro 

Ok.

Your model can be summarized as 

diff(G[dp](T), T) = psi[1]*F[d](T)+psi[2]*L[d](T)+psi[3]*E[c](T)-delta[3]*G[dp](T)
=
diff(G[dp](T), T) = F(T)-delta[3]*G[dp](T);

where
F(T) = psi[1]*F[d](T)+psi[2]*L[d](T)+psi[3]*E[c](T)

You say that F[d], L[d] and E[c] are indeed functions of time.
Let's assume that a correct fitting is such that psi[1], psi[2], psi[3], F[d](T), L[d](T) and E[c](T) get these values

psi[1] = p[1], psi[2] = p[2], psi[3] = p[3], F[d](T)=A(T), L[d](T)=B(T), E[c](T)=C(T)

Thus this solution is such that 

p[1]*A(T) + p[2]*B(T) + p[3]*C(T)

and any 6-uple 

(psi[1], psi[2], psi[3], F[d](T), L[d](T), E[c](T) )

which verifies 

psi[1]*F[d](T)+psi[2]*L[d](T)+psi[3]*E[c](T) = p[1]*A(T) + p[2]*B(T) + p[3]*C(T)

gives an exacyly as good fit than 

psi[1] = p[1], psi[2] = p[2], psi[3] = p[3], F[d](T)=A(T), L[d](T)=B(T), E[c](T)=C(T)

Thus your parameters are not identifiable.

If F[d], L[d] and E[c] arr known functions of time the problem remains the same: if F[d], L[d] and E[c] are not linearly independent there is an infinity of triples (psi[1], psi[2], psi[3]) which give exactly fit, meaning the psi are not identifiable.
Can you give the expressions of F[d], L[d] and E[c]? Or at least there values at the different years?

More of this you said nothing about the initial condition on G[dp](T)

Please enlighten us

• What does (for instance) year 2022 = 3, 023, 744. 14 mean?
  Are 3 quantities collected each year... or only one (GDP) and do you use commas to separate thousands?
• So, what are the data?
• What do you mean by "Fit the following model to the data" ?
  Do you mean "solve the ode

  diff(G[dp](T), T) = psi[1]*F[d](T)+psi[2]*L[d](T)+psi[3]*E[c](T)-delta[3]*G[dp](T);

  for some initial condition
  • Which one?  
    G[dp](some_unknown_date) = 1023744 ?
  • What about the parameters 

    delta[3] = .19, psi[1] = .63, psi[2] = .84, psi[3] = .72
    

    Are these the values to use in the previous ode?
    Or are there only indicative values (initial values) of these four parameters you aimed to assess through some minimization procedure?
  • Are F[d], L[d] andE[c] functions of time or constants?
     

Perhaps you think that the people you're asking for help just have to muddle through your draft, and that any effort on your part is superfluous?

@vv 

Your last remark is interesting.
There are no characters with more than 2 bytes in French (more precisely I didn't find any), that is why I asked @erik10 to check my procedure.

PS: I've just browsed the link you provided and voted for your answer.

... that you try to solve a system of 7 equations in only 4 unknowns?