@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

An advice for your future questions; avoid inserting the content of your uploaded file when it is that lond: insert only the link.

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

to your last question  here
 

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:
 

WHOEVER YOU ARE, PLEASE STOP CONVERTING THIS COMMENT INTO AN ANSWER.
THANK YOU


What is sigma_t[j]?
As you wrote it sigma_t[j] is a name which doesn't depend on any quantity, as you can verify doing this

indets(sigma_t[1], name)
                               {}

So:

diff(sigma_t[1], p[1]);
diff(sigma_t[1], E_0[90]);
diff(sigma_t[1], epsilon_dot);
                               0
                               0
                               0

Definining sigma_t does not imply defining sigma_t[j].

Finally:

1. Define sigma_t[j] correctly (or make sigma_t to depend on on some index j).
    

2. And once done define g should be defined this way

   g:= add(( sigma[j]-sigma_t[j])^2,j=1..10):
   

Infering what you want to achieve...

(in relation to your previous question), it seems to me that your problem could be stated this way:

•  sigma[j] is the outcome of an experiment done on some (viscoelastic?) material submitted to a loading epsilon.
   
• You have a physical model sigma_t  : true_strain ---> sigma parameterized by a vector p =(p[1], p[2], p[3]) which is aimed to explain the results observed duting 10 experiments whose input-output couples are ( true_strain[j], true_strain[j] ), j=1..10.
   
• Then g appears as the residual sum of squares between the observed outcomes sigma[j] and their predicted values sigma_t( true_strain[j] ).
   
• Would your goal be to find the vector p which minimizes g?
  If it is so the attached file will show you how to set correctly your problem.
   

  >  restart: >  sigma_t := epsilon -> E_0[90]*epsilon_dot*((epsilon/epsilon_dot)-(add(p[i], i = 1 .. 3)*(epsilon/epsilon_dot))+(add(p[i]*tau[i],i=1..3))-(add(p[i]*tau[i]*exp(-(epsilon/(epsilon_dot*tau[i]))),i=1..3))) (1) >  True_strain := Vector(10, symbol=epsilon): Sigma := Vector(10, symbol=sigma): >  g := add( (Sigma - sigma_t~(True_strain))^~2 ):   >  data   := indets(Sigma) union indets(True_strain); params := indets(g, name) minus data minus {seq(p[i], i=1..3)} (2) >  # Once data are known numerically, you must give "params" numeric values # # Then use something like Optimization:-NLPSolve(g) to find p Once data are known numerically, you must give "params" numeric values Then use something like Optimization:-NLPSolve(g) to find p Here is a notional example >  num_params := params =~ 1 (3) >  Observed_strain := Vector(10, i -> i); # Hypothesis : True_p := {seq(p[i]=i, i=1..3)}: # Then noisy values of observed sigma could be Observed_Sigma  := Vector(10, i -> eval( sigma_t(Observed_strain[i]), num_params union True_p ) + rand(-1. .. 1.)()); (4) >  g := add( (Observed_Sigma - eval(sigma_t~(Observed_strain), num_params))^~2 ); (5) >  Optimization:-NLPSolve(g, initialpoint=True_p, assume=nonnegative) (6) >

 

Download A_least_square_problem_maybe.mw

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   /     

Download Fibonacci.mw

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

@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

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

Download Boxplot_of_weighted_data.mw



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

 

Download Boxplot_of_real_weighted_data.mw

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

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

Download display_formula_mmcdara.mw

(UPDATED: please do not turn this comment into an answer)


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 doesn't look to the image you sent before:
 

Download solA_mmcdara.mw


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):
convert(add(s) + J^2*(-add(r)), horner, 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] /     /

@Andiguys 

# 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  

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

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

Point 1: Your differential system contains 8 parameters you have to fix before solving it numerically 

ind := convert(indets(ODE union IC, name) minus {eta}, list);
         [Le, Pr, alpha, N[b], N[t], a[1], a[2], a[3]]

Point 2: When you write

dsolve(eval(ODE union IC, a[i]=a), ...

what index "i" are you refering to?

Point 3: the output of your procedure finda is 

eval(f(eta),g(eta),h(eta),tt(10));

But eval needs two and only two argiments: what does this command of yours is meant to do?

Point 4: if you meant this

eval([f(eta),g(eta),h(eta)],tt(10));

the output of finda is a list with 3 elements, so 

fsolve(finda(a)-1, a=0..1);

is a nonsense because a list cannot be equal to a scalar ("1").

Look here for more explanations Shoot_Blasius-mmcdara.mw

See help(LinearAlgebra[Eigenvectors]): the output of Eigenvectors is a sequence val, vec wher val is a vector containing the eigenvalues and vec the matrix of eigenvectors.
To simplify the output of Eigenvectors do this

simplify([Eigenvectors(L)])

# or

val, vec := Eigenvectors(L):
simplify(vec)

Here is a more synthetic expression of the eigenvalues and eigenvectors:
synthetic.mw