7 years, 4 days

## @NanaP  I'm not sure I underst...

@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...

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

## You're right...

@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 thes...

@acer

Thank you acer for all these informations.

## Seems you ask two different question, do...

@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:

The Beta Subjective Distribution

 (1)

 (2)

 (3)

 (4)

 (5)

 (6)

 (7)

 > plot(f, x=Min..Max, gridlines=true)
 > mom := r -> evalf(Int(x^r*f, x=Min..Max)); _Mean     := mom(1); _Variance := mom(2)-_Mean^2; _Skewness := mom(3) / mom(2)^(3/2); _Kurtosis := mom(4) / mom(2)^2: Optimization:-Maximize(f, x=Min..Max); _Mode := eval(x, %[2]); Q := z -> evalf(Int(f, x=Min..z)): p := 0.25; _Q__||p := fsolve('Q(z)'=p, z=`if`(p < 0.5, Min, Max));
 (8)
 >

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.

The Beta Subjective Distribution

 (1)

 (2)

 (3)

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).

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-...

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.

## You're righ;...

@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 exp...

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

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

 > # A maple expression PutMeInTheBox = `#mrow(msubsup(mo("Σ"),mrow(mi("k"),mo("=120")),mo("1300")),mo("i",mathcolor="white"),msup(mrow(mo("("),mo("2"),mo("i",mathcolor="white"),mi("x-k"),mo(")")),mi("k")))`
 (1)
 >
 >

## Can't you provide a Maple code for u...

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 m...

@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:

 >
 >
 (1)
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 (2)
 >
 (3)
 >
 (4)
 >
 (5)
 >
 (6)
 >
 >
 > W     := unapply(W, Y);Theta := unapply(Theta, Y); U     := unapply(U, Y);Phi   := unapply(Phi, Y); F     := unapply(F, Y);T     := unapply(T, Y); G     := unapply(G, Y);H     := unapply(H, Y);
 (7)
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >

 >
 > Zs :=  [z1, z2, z3, z4, z5, z6, z7, z8, z9, z10, z11, z12, z13, z14,z15, z16,z17,z18,z19,z20,z21, z22, z23, z24,z25, z26]: print~(Zs):
 (8)
 > IZs := indets(Zs) minus {M}; numelems(Zs), numelems(IZs);
 (9)
 > Zsol := proc(Ma) fsolve(eval(Zs, M=Ma)) end proc: # example Zsol(1)
 (10)
 > F := proc(Y, Ma) eval(sum(b[i]*Y^i,i=0..7), Zsol(Ma)) end proc: T := proc(Y, Ma) eval(sum(c[i]*Y^i,i=0..7), Zsol(Ma)) end proc: G := proc(Y, Ma) eval(sum(d[i]*Y^i,i=0..4), Zsol(Ma)) end proc: H := proc(Y, Ma) eval(sum(h[i]*Y^i,i=0..4), Zsol(Ma)) end proc:
 >
 > plot([seq(T(Y, Ma), Ma = 1 .. 3)], Y = 0 .. 1, axes = boxed, color = [red, green, blue], legend = `~`[cat]("M=", [`\$`(1 .. 3)]), thickness = 2, labels = [Y, T])
 >
 >
 >
 (11)
 >
 (12)

 > Cf := 'Cf': Nu := 'Nu': plot(   [seq(a1*(1+1/bet)*diff(F(Y, Ma), Y), Ma = 1 .. 3)]   , Y=0..1   , color=[red, green, blue]   , legend = `~`[cat]("M=", [`\$`(1 .. 3)])   , labels=[Y, "Cf"] ); print(): plot(   [seq((-a4-Rd)*diff(T(Y, Ma), Y), Ma = 1 .. 3)]   , Y=0..1   , color=[red, green, blue]   , legend = `~`[cat]("M=", [`\$`(1 .. 3)])   , labels=[Y, "Nu"] )

 >
 >
 >
 >

## @ogunmiloro Ok.Your model can be su...

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)

## Could you make some presentation effort?...

• 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 ?
```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?