8 years, 43 days

## Two ways (at least)...

The first one consists in "rescaling the abscissas" of a classical ErrorPlot; the second one (which I prefer for it is more versatile) is to create "from scratch" your own error plot procedure.

ErrorPlot.mw

Example of what the second approach can produce:

As a rule I prefer  to define my own visualization procedure (of statistical data) instead of using those provided by the Statistics package.
This doesn't mean these latter are not "good" but only that they are not versatile enough to fit my requirements

## Any constant function u is a solution...

Let us consider now this much simpler problem

```(a*x+b*y)*diff(u(x, y), x) + (c*x+d*y)*diff(u(x, y), y) = 0
```

For this problem to have a solution you must have a=-d, which gives

```                         1    2           1    2
u(x, y) = - b y  + a x y - - c x
2                2
```

So your problem could have a non trivial solution (aka a constant one) only for some values of values of alpha, b, c, d.
Here is an example of a non trivial solution for the case alpha=1, c=0, d=0 compatibility_2.mw

Another point of view
Your equation writes F(x, y, z, w) ° gradient(u(x, y, z, w) where F is field of directions in space [x, y, z, w] and ° denotes the dot product.
This direction field defines a family of lines such that the tangent to each of them is confounded with the direction of the field at the contact point.
Defining x(s), ..., w(s) as the cosine of this family gives:

```                 d
--- x(s) = alpha (b w(s) + y(s))
ds
d
--- y(s) = -b w(s) + x(s) + z(s)
ds
d
--- z(s) = -c y(s)
ds
d
--- w(s) = d (y(s) - x(s))
ds
```

Integrating this system gives the expression of the characteristic lines.
Then u is an hypersurface made of these characteristic lines.

Look here  (characteristic system can be solved numerically but a numerical solution is provided as an illustration):

 > restart
 > E0 := alpha*(y+b*(w))*diff(u(x,y,z,w),x)+(x+z-b*(w))*diff(u(x,y,z,w),y)-c*y*diff(u(x,y,z,w),z)+d*(y-x)*diff(u(x,y,z,w),w)=0;
 (1)
 > eo := [op(lhs(E0))]
 (2)
 > ed := map2(op, -1, eo)
 (3)
 > ec := eo /~ ed;
 (4)
 > ew := map2(op, -1, ed)
 (5)
 > # Classical mathematical representation zip((s1, s2) -> cat(d, s1)/s2 = ds, ew, ec)
 (6)
 > # Maple translation zip((s1, s2) -> diff(s1(s), s) = s2, ew, ec); sys := lhs~(%) =~ eval(rhs~(%), ew =~ ew(s)): print~(sys):
 (7)
 > Characteristic_lines := dsolve(sys, ew(s)): whattype(%); numelems(%%);
 (8)
 > # Numerical appoximation params := convert(indets(E0, name) minus convert(ew, set), list): data := params =~ [seq(rand(0. .. 1.)(), i=1..numelems(params))];
 (9)
 > numsys := eval({sys[], (ew(0) =~ [seq(rand(-1. .. 1.)(), i=1..4)])[]}, data)
 (10)
 > numsol := dsolve(numsys, numeric)
 (11)
 > with(plots): col := [red, magenta, cyan, blue]: display(   seq(     plots:-odeplot(numsol, [s, ew[i](s)], s=0..4, color=col[i], legend=ew[i](s))     , i=1..4   ) )

A more general approach

 > AllParams := [params[],  S, seq(C__||i, i in ew) ]; numsys    := {sys[], seq(ew[i](S) = C__||(ew[i]), i=1..4)}
 (12)
 > numsol := dsolve(numsys, numeric, parameters=AllParams)
 (13)
 > # As a first illustration only the initial value S is changed while all the # other parameters are ketp constant (same value 1, but you can change them) V := proc(sigma)   numsol(parameters=[1\$4, sigma, 1\$4]):   display(     seq(       plots:-odeplot(numsol, [s, ew[i](s)], s=0..4, color=col[i], legend=ew[i](s))       , i=1..4     )   ) end proc: # Explore(V(sigma), parameters=[sigma=0. .. 10.])
 > # As a first illustration only the initial value S is changed while all the # other parameters are ketp constant (same value 1, but you can change them) V := proc(sigma)   numsol(parameters=[1\$4, sigma, 1\$4]):   odeplot(numsol, [s, x(s)], s=0..4, color=blue, legend=x(s)) end proc: #display( seq(V(sigma), sigma=0..5, 0.1), view=[default, -2..2])
 > V := proc(sigma)   numsol(parameters=[1\$4, sigma, 1\$4]):   odeplot(numsol, [s, x(s), y(s)], s=0..4, color=blue) end proc: #display( seq(V(sigma), sigma=0..5, 0.1))
 > V := proc(sigma)   numsol(parameters=[1\$4, sigma, 1\$4]):   odeplot(numsol, [x(s), y(s), z(s)], s=0..4, color=blue) end proc: #display( seq(V(sigma), sigma=0..5, 0.1))

## @Aixleft aimbot  Here is a more re...

@Aixleft aimbot

Here is a more readable code where some errors have been corrected.

Nevertheless it remains 4 expressions which "come from nowhere" and which I wasn't capable to check if the are correct or not: N11, N22, N33, N44 (it seems that you rewote yourself the expressions of N1, N2, N3, N4... Am I right?)

If it is so here are the expressions you are looking for:

```
2    1   /   (1/2)    \  (1/2)  3    / 4\
tr :=    2 - v  + ---- \3 5      + 5/ 5      v  + O\v /
1200

1   /   (1/2)    \  (1/2)  3    / 4\
det :=       1 - ---- \7 5      + 5/ 5      v  + O\v /
1200
```

## mtaylor...

@felixiao

The output of taylor(...) is algebraic_expression+O(..) : the presence of O(..) causes the error you get.
The simplest thing to do is to use mtaylor instead whose output is an algebraic_expression.

A lengthy way is to write convert(taylor(...), polynom) instead of taylor(...).

MMSE_4COS_ok.mw

## try & catch...

Illustrative example

```restart:

PLSQ := (x, y) -> 1/x+1/y;

for i from -1 to 1 do
for j from -1 to 1 do
try
print(i, j, PLSQ(i, j)):
catch:
print(i, j, "division by 0 not allowed")
end try
end do:
end do:
-1, -1, -2
-1, 0, "division by 0 not allowed"
-1, 1, 0
0, -1, "division by 0 not allowed"
0, 0, "division by 0 not allowed"
0, 1, "division by 0 not allowed"
1, -1, 0
1, 0, "division by 0 not allowed"
1, 1, 2
```

## See my comment...

To know what are the methods Optimization uses, just do this

```help("Methods Used by the Optimization Package")
```

More specifically, concerning Statistics:-NonlinearFit run this command

```help("Regression Options")
```

From this help page:

## First and second proofs...

The proof of the first formula is very simple but I rewrite F into F2 forreasons that will appear clearly during the the proof of your second formula.

Formula 1:

```F := unapply(rsolve({F(1) = 1, F(2) = 1, F(n + 1) = F(n) + F(n - 1)}, F(n)), n):

rel2use := Phi = op([1, -1, 1], F(n)):
rel2use := { rel2use, C = op(1, F(n)) / rhs(rel2use)^n }

/    1  (1/2)        1  (1/2)   1\
{ C = - 5     , Phi = - 5      + - }
\    5               2          2/

# FORMULA 1

F2 := n -> C * (Phi^n -(1-Phi)^n):

ToProve_1 := F2(n + 1)^2 - F2(n)*F2(n+2) = (-1)^n:

expand~(ToProve_1):
simplify(eval(%, rel2use))
n       n
(-1)  = (-1)
```

Formula 2:

```# Preliminary:

G2 := n -> arccot(``(F2)(n));
G2 := n -> arccot((F2)(n))

ToProve_2 := G2(2*n) - (G2(2*n + 1) + G2(2*n + 2)) = 0;
arccot((F2)(2 n)) - arccot((F2)(2 n + 1)) - arccot((F2)(2 n + 2)) = 0

# From recurrence formula, replacing n+1 ny 2*n+2 and F by F2:

recur := ``(F2)(2*n+2) = ``(F2)(2*n+1) + ``(F2)(2*n);
(F2)(2 n + 2) = (F2)(2 n + 1) + (F2)(2 n)

# One gets:

subs(recur, ToProve_2);
arccot((F2)(2 n)) - arccot((F2)(2 n + 1)) - arccot((F2)(2 n + 1) + (F2)(2 n)) = 0

# Which has symbolic form:
# (U > 0, V > 0, U < V)

SymbolicForm := arccot(U) - arccot(V) - arccot(U+V)
arccot(U) - arccot(V) - arccot(U + V)

# Proof

G2 := n -> arccot(F2(n)):

recur := F2(2*n+2) = F2(2*n+1) + F2(2*n):

ToProve_2 := simplify(G2(2*n) - (G2(2*n + 1) + G2(2*n + 2)) = 0):

SymbolicForm := -(1/2)*Pi + arctan(U) - arctan(V) - arctan(U+V)
1
- - Pi + arctan(U) - arctan(V) - arctan(U + V)
2

# Thus

SymbolicForm := -(1/2)*Pi + arctan((U-V)/(1+U*V)) - arctan(U+V):

# and next

SymbolicForm := simplify( -(1/2)*Pi + arctan(((U-V)/(1+U*V) - (U+V)) / (1+(U-V)/(1+U*V)*(U+V))) )
/  / 2          \ \
1            |V \U  + U V + 2/ |
- - Pi - arctan|-----------------|
2            | 2          2    |
\U  + U V - V  + 1/

FocusOn := denom( op(1, indets(SymbolicForm, function)[]) )
2          2
U  + U V - V  + 1

FocusOn := eval(FocusOn, {U=F2(2*n), V=F2(2*n+1)}):
FocusOn := eval(FocusOn, rel2use):
FocusOn := simplify(FocusOn) assuming n::integer;
(4 n)                     (4 n + 1)
1  /  1  (1/2)   1\        1 /  1  (1/2)   1\
- -- |- - 5      + -|      + - |- - 5      + -|
10 \  2          2/        5 \  2          2/

(4 n)
1  /  1  (1/2)   1\       (1/2)
+ -- |- - 5      + -|      5
10 \  2          2/

FocusOn := ``(Y^(4*n)) * expand(eval(FocusOn, (-(1/2)*sqrt(5)+1/2)=Y) / Y^(4*n))
/ (4 n)\ /  1    1     1   (1/2)\
\Y     / |- -- + - Y + -- 5     |
\  10   5     10       /

eval(FocusOn, Y=-(1/2)*sqrt(5)+1/2);
0
# Then

'SymbolicForm' = op(1, SymbolicForm) + arctan(something/``(0))
1            /something\
SymbolicForm = - - Pi + arctan|---------|
2            \   (0)   /

something := numer( op(1, indets(SymbolicForm, function)[]) )
/ 2          \
V \U  + U V + 2/

# Obviously 'something' is strictly positive, thus
#
# arctan(something/``(0)) = Pi/2;
#
# and thus

'ToProve_2' = ``(-Pi/2) + ``(Pi/2)

/  1   \    / 1   \
ToProve_2 = |- - Pi| + |  - Pi|
\  2   /    \ 2   /

```

PS I believe the second proof could be easier to catch if one uses the trigonometric representation of the golden number.

## @jalal " all calculators ...

@jalal

" all calculators are programmed to divide by n... it's a broad debate"
Not all the calculators are programmed to divide by n (look to Excel(*) which proposes to divide by n or by n-1).
And it is not a broad debate: the situation where you have to divide by n is extremely clear, as is the situation where you have to divide by n-1.
The only situation when you have to divide by n is when your sample represents the whole population. If it is not the case dividing by n lead to a biased estimator of the variance, while dividing by n-1 lead to an unbiased estimator.

The reason I used n-1 (without asking you which one of these two situations you were in) is a direct consequence of the way you construct data and Weights:

• As you randomly chose the elements of data from the range 4..20 this implicitely means that you draw a sample from a population whose indicuals are 4, 5, ..., 20. So data does not represents the whole population.
• The same reasoning applies to Weights which are randomly chosen.

Had you said "I have population made of W[1] individuals X[1], ..., W[N] individuals X[N]" , the situation would have been different and the division by n perfectly justified.

Excel(*) : Excerpt from  Variance Excel:

Excerpt from Population variance Excel

REMARK: Alternative ways to generate data are

```N := 10:

data := sort(combinat:-randperm([\$4..20])[1..N]);
data := sort(Statistics:-Shuffle([\$2..20])[1..N]);

CodeTools:-Usage( sort(Statistics:-Shuffle([\$2..20])[1..N]), iterations=10^5):
memory used=4.73KiB, alloc change=0 bytes, cpu time=46.96us, real time=47.03us, gc time=2.25us

CodeTools:-Usage( sort(combinat:-randperm([\$4..20])[1..N]), iterations=10^5):
memory used=2.48KiB, alloc change=0 bytes, cpu time=15.79us, real time=15.80us, gc time=1.14us

```

## Do you mean......

@jalal

...that you have weighted data and that you want to draw a boxplot of these data?

if it is so the trick is to build an unweighted sample from the data and their weights and draw the corresponding boxplot.
This is likely a hammer to kill a fly, but the method which consists in coding your own bowplot is far more complex (see Add-on 2 in attached file number two)

If the weights are strictly positive integers (which is the simpler case) use this:

 > restart
 > with(Statistics):
 > N       := 10: data    := Sample(Normal(0, 1), N):
 > # non-weighted data evalf[5](FivePointSummary(data));
 (1)
 > # weighted data weights := [ seq(rand(1..10)(), n=1..N) ]; sample := `<,>`( seq(data[n]\$weights[n], n=1..N) ); evalf[5](FivePointSummary(sample));
 (2)
 > BoxPlot([data, sample], axis[1]=[tickmarks=[1="unweighted", 2="weighted"]])
 >

If weights are strictly positive real numbers (a little bit more complex case) use this file
(look more specifically  the Maple issue of Add-on 3).

 > restart
 > with(Statistics):
 > Seed := 171803456472246; randomize(Seed);
 (1)
 > N       := 10: data    := sort(Sample(Normal(0, 1), N)): weights := [ seq(rand(0. .. 10. )(), n=1..N) ]:
 > # non-weighted data evalf[5](FivePointSummary(data)); print(): interface(rtablesize=N+1): VW := evalf[5](`<,>`(`<|>`("Values", "True weights"), `<|>`(data^+, < weights >)));
 (2)

FROM WEIGHTED TO UNWEIGHTED DATA

 > prob := weights /~ add(weights): PT   := RandomVariable(ProbabilityTable(prob)): UseHardwareFloats := false: # The larger "LargeNumber", the closer the weights in sample to the real weights "weights" LargeNumber := 10^5: DataNumbers := round~(Sample(PT, LargeNumber)): sample := Vector(LargeNumber, i -> data[DataNumbers[i]]): evalf[5](FivePointSummary(sample)); print(): TS := Tally(sample): ST := sort(evalf[4](lhs~(TS)) =~ rhs~(TS), key = (x -> lhs(x))): # What are the weights we generate? SW := evalf[5]( rhs~(ST) /~ (LargeNumber/add(weights)) ): `<|>`(VW, Vector(["Weights in sample", SW[]]))
 (3)
 > BoxPlot([data, sample], axis[1]=[tickmarks=[1="unweighted", 2="weighted"]])

CONSTRUCTION OF AN UNWEIGHTED SAMPLE

IMO more complex than what I did above, but could be maybe an option when the size of data is large

 > Unweight := proc(Data, Weights, LargeNumber)   local prob, PT, DataNumbers, sample:       UseHardwareFloats := false:   prob        := weights /~ add(weights):   PT          := RandomVariable(ProbabilityTable(prob)):   DataNumbers := round~(Sample(PT, LargeNumber)):   sample      := Vector(LargeNumber, i -> data[DataNumbers[i]]):   return sample end proc:
 > BoxPlot([data, Unweight(data, weights, 10^5)], axis[1]=[tickmarks=[1="unweighted", 2="weighted"]])

BOXPLOT RECONSTRUCTION OF WEIGHTED DATA

IMO more complex than what I did above, but could be maybe an option when the size of data is large

 > WBoxPlot := proc(Data, Weights, MyColor)   local W_mean, W_median, W_Q25, W_Q75, W_deciles,         W_upper, W_lower, W_outlier, W_boxplot,         cumprob, k:   uses plots, plottools:   W_mean    := Mean    (Data, 'weights'=Weights):   W_median  := Median  (Data, 'weights'=Weights):   W_Q25     := Quantile(Data,  0.25, 'weights'=Weights):   W_Q75     := Quantile(Data,  0.75, 'weights'=Weights):   # This is the natural command you might thinkto to compute deciles.   # But as shown in Add-on 3 it returns wrong results.   # W_deciles := seq(Quantile(Data,  d, 'weights'=Weights), d=0.1..0.9, 0.1):   #   # So I use instead   cumprob   := CumulativeSum(Weights) / add(Weights):   W_deciles := NULL:   for k from 0.1 to 0.9 by 0.1 do     W_deciles := W_deciles, evalf[5]( data[ListTools:-BinaryPlace(cumprob, k)+1] )   end do;   W_upper   := min(max(Data), W_Q75+3/2*(W_Q75-W_Q25)):   W_lower   := max(min(Data), W_Q25-3/2*(W_Q75-W_Q25)):   W_outlier := Select((t -> is(t < Q_Lower or t > W_upper)), Data):   W_boxplot := display(     rectangle([-1, W_Q25], [1, W_Q75], color=MyColor)     , plot([[-1, W_median], [1, W_median]], color=black)     , plot([[0, W_Q75], [0, W_upper]], color=black)     , plot([[0, W_Q25], [0, W_lower]], color=black)     , plot([[-0.5, W_upper], [0.5, W_upper]], color=black)     , plot([[-0.5, W_lower], [0.5, W_lower]], color=black)     , pointplot([seq([0, W_deciles[i]], i=1..9)], symbol=circle, symbolsize=15, color=black)   ):   return W_boxplot end proc:
 > use plots, plottools in   display(     BoxPlot([data, sample])     , translate(scale(WBoxPlot(data, weights, "ForestGreen"), 0.75/2, 1), 3, 0)     , axis[1]=[tickmarks=[1="unweighted", 2="weighted", 3="weighted like boxplot"]]   ) end use

A MAPLE ISSUE

Maple uses several methods to compute quantile (referenced by an integer from 1 to 8, see the corresponding help page).
The default method has number 7 and gives unexpected results when unequal weights are used.
This issue affects Percentile and Decile too.

So be extremely carefull while using these procedures with unequal weights

 > printf("This value is correct\n"): Median    (data, 'weights'=weights); printf("\nThese three values (method=default) are wrong\n"): Quantile  (data,  0.5, 'weights'=weights); Percentile(data,  50, 'weights'=weights); Decile    (data,  5, 'weights'=weights); # This is the reason why percentiles are not correctly placed on the green boxplot
 This value is correct
 These three values (method=default) are wrong
 (4)
 > printf("\nThese three values (method=1) are correct\n"): Quantile  (data,  0.5, 'weights'=weights, 'method'=1); Percentile(data,  50, 'weights'=weights, 'method'=1); Decile    (data,  5, 'weights'=weights, 'method'=1); print(): printf("\nThese three values (method=2) are correct\n"): Quantile  (data,  0.5, 'weights'=weights, 'method'=2); Percentile(data,  50, 'weights'=weights, 'method'=2); Decile    (data,  5, 'weights'=weights, 'method'=2);
 These three values (method=1) are correct
 These three values (method=2) are correct
 (5)
 > # The computations are correct for the unweighted sample: Median    (sample); Quantile  (sample, 0.5); Percentile(sample, 50); Decile    (sample, 5);
 (6)
 > # How to compute correctly deciles of weighted data cumprob  := CumulativeSum(prob): deciles := NULL: for k from 0.1 to 0.9 by 0.1 do   deciles := deciles, [k, evalf[5]( data[ListTools:-BinaryPlace(cumprob, k)+1] )] end do: `<,>`(< ':-Probability' | ':-Decile' >, convert([deciles], Matrix))
 (7)
 > # Check these values by using the unweighted sample deciles := NULL: for k from 1 to 9 do   deciles := deciles, [k, evalf[5]( Decile(sample, k) )] end do: `<,>`(< ':-Probability' | ':-Decile' >, convert([deciles], Matrix))
 (8)
 > # MAPLE BUG wrong_deciles := NULL: for k from 1 to 9 do   wrong_deciles := wrong_deciles, [k, evalf[5]( Decile(data, k, 'weights'=prob) )] end do: `<,>`(< ':-Probability' | ':-Decile' >, convert([wrong_deciles], Matrix))
 (9)
 > # Check that BoxPlot correctlycomputes deciles for unweihghted sample buw := BoxPlot(sample): op(1..-2, select(has, [op(buw)], POINTS)[2]): Deciles_from_BoxPlot := map(k -> [k, evalf[5](%[k][2])], [\$1..9]): `<,>`(< ':-Probability' | ':-Decile' >, convert(Deciles_from_BoxPlot, Matrix));
 (10)
 >

NOTE:
There exist several quantile estimators for weighted data and the one I used here is just the most commonly used and easier to understand.

See for instance
Akinshin

## If I'm not mistaken...

To moderators: please stop converting this comment into an answer, thank you

the expression you want to textplot is the result of the optimization of a Lenard-Jones model you use in your previous question?

if it is so the attached file contains a full reply (likely not the same data points the LJ model comes from).

 > restart;
 > pts := [[1.020408163, .9618194785], [1.021860259, .9591836735], [1.047329717, .9183673469], [1.077147904, .8775510204], [1.112217401, .8367346939], [1.153643677, .7959183673], [1.202786007, .7551020408], [1.224489796, .7394037725], [1.266731737, .7142857143], [1.350499968, .6734693878], [1.428571429, .6426647396], [1.463093584, .6326530612], [1.632653061, .5927592700], [1.638566951, .5918367347], [1.836734694, .5653094517], [2.024654644, .5510204082], [2.040816327, .5499100928], [2.244897959, .5415404104], [2.448979592, .5375577894], [2.653061224, .5364038532], [2.857142857, .5371174931], [3.061224490, .5390909888], [3.265306122, .5419306191], [3.469387755, .5453753581], [3.673469388, .5492483720], [3.759428515, .5510204082], [3.877551020, .5532080969], [4.081632653, .5571718264], [4.285714286, .5612391621], [4.489795918, .5653739340], [4.693877551, .5695506292], [4.897959184, .5737510616], [5.102040816, .5779621428], [5.306122449, .5821743709], [5.510204082, .5863807930], [5.714285714, .5905762856], [5.775983717, .5918367347], [5.918367347, .5943431266], [6.122448980, .5978993957], [6.326530612, .6014214453], [6.530612245, .6049104007], [6.734693878, .6083674439], [6.938775510, .6117937612], [7.142857143, .6151905094], [7.346938776, .6185587956], [7.551020408, .6218996660], [7.755102041, .6252141001], [7.959183673, .6285030098], [8.163265306, .6317672398], [8.219364040, .6326530612], [8.367346939, .6345736540], [8.571428571, .6371910423], [8.775510204, .6397830180], [8.979591837, .6423509843], [9.183673469, .6448962158], [9.387755102, .6474198713], [9.591836735, .6499230055], [9.795918367, .6524065794], [10.00000000, .6548714697]]:
 > model := a[1] + a[2]*((a[3]/Gamma)^a[4] - (a[5]/Gamma)^a[6]);
 (1)
 > K := unapply(model, Gamma); J := add((map2(op, 2, pts) -~ K~(map2(op, 1, pts)))^~2): opt := evalf[6]( Optimization:-NLPSolve(J, assume=nonnegative) );
 (2)
 > alias(STS = StringTools:-Substitute): string_model := convert(model, string); values       := evalf[3](opt[2])
 (3)
 > for n from 1 to 6 do   string_model := STS(string_model, convert(lhs(values[n]), string), convert(rhs(values[n]), string)) end do:
 > plots:-textplot([0, 0, cat("W__LJ = ", string_model)], axes=none,'font'=["helvetica","roman",25], size=[1000, 100])
 >

## There are different way to fix this...

In your procedure the same quantity m is computed each time the loop is executed. So when theloop ends m contains the last result (here 11).
If you want to procedure to return all the values of m, you must record them in some structure.

For instance (tst returns a list)

```tst := proc (M)
local i, m;
m := NULL;
for i from 2 to numelems(M)+1 do
m := m, M[1]-M[i-1]
end do;
return [m]
end proc:```

Or (tst returns a vector)

```tst := proc (M)
local i, m;
m := Vector(numelems(M)):
for i from 2 to numelems(M)+1 do
m[i-1] := M[1]-M[i-1]
end do;
return m
end proc:```

Or more simply;

`M[1] -~ M[2..-1]`

## @mz6687 Here is the solution but it...

@mz6687

Here is the solution but it doesn't look to the image you sent before:

 >
 >
 >
 >
 >
 >
 >
 > sol := solve({x2, x1, s2}, {m0, n0, lambda})
 (1)
 > valsol := simplify~({allvalues(sol)});
 (2)
 > nops(valsol);  # You gave two solutions, not one
 (3)
 > # Here is for instance solution 1 (solution 2 has the same form) print~(valsol[1]):
 (4)
 > # A simpler display Z := 'Z': R := select(type, indets(valsol), `^`)[] = Z: ToDisplay := collect~(subs(R, valsol[1]), Z): print~([(rhs=lhs)(R), ToDisplay[]]):
 (5)
 > ToDisplay
 (6)
 > # Or even simpler: ToDisplay := map(u -> lhs(u) = value(``(numer(rhs(u))) / denom(rhs(u))), ToDisplay): print~([(rhs=lhs)(R), ToDisplay[]]):
 (7)
 >

If you really want this "imaginary" form your imagedisplays (I don't find any interest in it), you can obtain it after some dymnastis:

```eval(lhs(R), op(1, lhs(R)) = -op(1, lhs(R))*J^2):
newR := op(2, (combine(simplify(%), radical) assuming J > 0)) = Y:
eval( collect(subs(R, valsol[2][1]), Z), Z=J*Y):

s, r := selectremove(has, [op(rhs(%))], J):
J * value(``(numer(op(1, %))) / denom(op(1, %))):

lambda = subs(J=I, %);
With = ( (rhs=lhs)(newR) )
1
lambda = - I (I a[1] + I a[2] + Y)
4
/                                             (1/2)\
|    /    2         2                       2\     |
With = \Y = \12 c  + 3 a[1]  - 6 a[1] a[2] + 3 a[2] /     /
```

## Obvious solution:...

```# First thing: construct TRC correctly

TRC_ok := eval(TRC, {alpha = (u -> 0.02*exp(0.4*u)), beta = (u -> 0.2e-1*(1-exp(.5*I2)))})
12
4.267000015 10   + 90 Q2 - 21.25 Q1

7
+ 4.000000000 10  I1 exp(0.4 I1)

+ 2000000000 I2 (0.02 - 0.02 exp(0.5 I2))

# Assuming Q2 is nonnegative, TRC_ok is obviously mimimum when Q2=0.

TRC_reduced := eval(TRC_ok, Q2=0);
12                            7
4.267000015 10   - 21.25 Q1 + 4.000000000 10  I1 exp(0.4 I1)

+ 2000000000 I2 (0.02 - 0.02 exp(0.5 I2))

# To be sure of what you are about to do verify what are the indeterminates of TRC_ok:

indets(TRC_ok, name)
{I1, I2, Q1, Q2}

# As you commented C1, C2 and C3, the Minimize command cannot work

C1 := 4e8 >= Q1,  Q1 >= .6e7;
C3 := I1 <= I2;
8      6
Q1 <= 4 10 , 6 10  <= Q1
I1 <= I2
```

Next:

```# Do you really want to write this:
#  C2 := if Q1= 4*(10)^(8) then Q2 =0 else Q2 = (x*d*h)-(delta*Q1*h)
# which meand that Q2 is always equal to (x*d-delta*Q1)*h excepted when Q2=4e8 where Q1=0 ???
#
# As you plane to use a numerical algorithm (Minimize) the probability that Q2=4e8 at some
# iteration is exactly null... unless the minimum verifies Q1=4e-8 (limit of the first C1
# condition).
```

Can we find "by hand" the minimum you are looking for?
TRC_reduced is obviously minimal when Q1 is maximal (Q1=4e8).

```# The expression of TRC_reduced shows it is minimal when Q1 is maximal, thus

evalf(eval(TRC_reduced, Q1=4e8));
12         7                        9
4.258 10   + 4.0 10  I1 exp(0.4 I1) + 2.0 10  I2 (0.02 - 0.02 exp(0.5 I2))
```

As a simple plot shows, the derivative wrt I2 of this expression is negative for 0 <= I2 <= 6.
Thus the minimum of TRC_reduced corresponds to Q1=4e8 and  I2=6.
It remains

```evalf(eval(TRC_reduced, [Q1=4e8, I2=6]));
12         7
4.254 10   + 4.0 10  I1 exp(0.4 I1)
```

... which is obviously minimum for I1=0 (a value which verifies I1 <= I2):

```evalf(eval(TRC_reduced, [Q1=4e8, I2=6, I1=0]));
12
4.254 10

```

## A bins-and-counts histogram...

(By the way I think you already asked  this same question and I already answered it. As I don't have time to verify this, please forgive me I'm wrong).

Whatever, here is a way to build an histogram from couples (bins, counts).
I assume your data are a list of equalities bin=count, where bin=a..b is a range and count is an integer.
For other structures, for instance a two-column matrix whose column 1 contains the bins and column 2 the counts, proceed as explained at the end oftheattached file.

Bins_and_counts_histogram.mw

I remain at your disposal for further help.

Remark: in the arrached file STEP 2 shows how to get the bins and counts of raw data while using TallyInto.

## Why are you talking about regression?...

If you want to estimate the size and shape parameters of the Weibull random variable the data S could bea sample drawn from, you just have to do this

```mle := MaximumLikelihoodEstimate(Weibull(a, b), S);
mle := [a = 129.944334168688, b = .9301095062]
```

If you want to build the corresponding random variable for any other purpose, I advice you to do this

```W   := RandomVariable(Weibull(a, b)):
mle := MaximumLikelihoodEstimate(W, S):
W__opt := Specialize(W, mle):

# Example of use

N := numelems(S);
QuantilePlot(Sample(W__opt, N), S):
```
 1 2 3 4 5 6 7 Last Page 1 of 55
﻿