mmcdara

635 Reputation

10 Badges

3 years, 41 days

MaplePrimes Activity


These are answers submitted by mmcdara

The problem is such general, with such many potential situations, that no universal answer can be given (it's my opinion).
So, if you have a particular example in mind, you should load the mw file by using the big green arrow in order a more suited answer might be delivered.

Nevertheless here are two examples of what it is possible to do


ATTENTION : I am not sure that 'allvalues' always returns the solutions in the same order order... and this is essential here to ensure that the curves are not tangled.
This is a point I will verify and, possibly modify the code accordingly

 

restart:

interface(rtablesize=10):

# lets's start with a simpple case
# Suppose that x and y are two functions of z of KNOWN explicit forms
# Example

x := z -> cos(z):
y := z -> z*sin(z):

# then the "solution curve" C(z) = { (x(z), y(z)), z=a..b } can be drawn with

plot([x(z), y(z), z=0..20*Pi]);

 

# Now a harder case:
# x and y are two functions of z of UNKNOWN explicit forms
# Example

vars := [x, y, z, t];
randpoly(vars, degree = 2);

[x, y, z, t]

 

-73*t^2-56*x^2-62*y^2+97*z^2+87*x

(1)

sys := [ seq(randpoly(vars, degree = 2)=0, n=1..3) ]:
print~(sys):

71*x*z-17*y^2-75*y*z-44*t+80*x-82 = 0

 

74*t*y+37*t*z-92*y^2+6*y*z+72*z^2+75*x = 0

 

-47*t*z-29*x^2+95*x*y+11*y*z-49*z^2-8 = 0

(2)

# Exact solutions

ExactSolutions := solve(evalf(sys),vars):

# The solution is made of 1 functional solution and 10 "numeric" 4-uples

NumberOfSolutions := numelems(ExactSolutions);

for n from 1 to NumberOfSolutions do
 # printf("List of functions the solution %2d contains %a\n", n, select~(has, rhs~(ExactSolutions[n]), 'function') )
  printf("List of functions the solution %2d contains %a\n", n, numelems~(indets~(rhs~(ExactSolutions[n]), function)) )
end do:

11

 

List of functions the solution  1 contains [1, 1, 0, 1]
List of functions the solution  2 contains [0, 0, 0, 0]
List of functions the solution  3 contains [0, 0, 0, 0]
List of functions the solution  4 contains [0, 0, 0, 0]
List of functions the solution  5 contains [0, 0, 0, 0]
List of functions the solution  6 contains [0, 0, 0, 0]
List of functions the solution  7 contains [0, 0, 0, 0]
List of functions the solution  8 contains [0, 0, 0, 0]
List of functions the solution  9 contains [0, 0, 0, 0]
List of functions the solution 10 contains [0, 0, 0, 0]
List of functions the solution 11 contains [0, 0, 0, 0]

 

# Only solution 1 is functional with x, y and t expressed as functions of z:

ExactSolutions[1][3]

z = z

(3)

# Here are numerical values of the couple (x, y)  if z=0 (for instance)

evalf~(subs~(z=0, rhs~(ExactSolutions[1][1..2])))

[[.4652284527, .3230260837]]

(4)

# doing so for fifferent values of z enables ploting the "solution curve" C(z) = { (x(z), y(z)), z=a..b }
# for some real range a..b
# Example

f := u -> evalf~(subs~(z=u, rhs~(ExactSolutions[1][1..2]))):

Cxy := [seq(f(u), u in seq(0..10, 0.1))]:

plot(Cxy, labels=['x(z)', 'y(z)']);



Cxyz := [seq([f(u)[], u], u in seq(-10..10, 0.1))]:

PLOT3D(CURVES(Cxyz), AXESLABELS('x(z)', 'y(z)', 'z'))

 

 

# The problem is that x(z) is not unique and that among all the solutions some are complex

[ evalf(allvalues(subs(z=0, rhs(ExactSolutions[1][1])))) ];
numelems(%);

# and neither is y(z)
[ evalf(allvalues(subs(z=0, rhs(ExactSolutions[1][2])))) ];
numelems(%);

[.4652284527, 40.93176235, -0.2722500350e-1-0.2586203384e-1*I, -.6917623044+.5507283075*I, -.6917623044-.5507283075*I, -0.2722500350e-1+0.2586203384e-1*I]

 

6

 

[.3230260837, 12.49701637, -1.634238234+1.536633637*I, -.2856782808+.1087988986*I, -.2856782808-.1087988986*I, -1.634238234-1.536633637*I]

 

6

(5)

# Test if Cxyz(s) is real before doing the plot.

realCxyz := proc(a, b, step, s, verbose)
   local g, c, u, h:

   g := (u, s) -> evalf~(allvalues(subs~(z=u, rhs~(ExactSolutions[1][1..2])))[s]):

   c := NULL:
   for u from a to b by step do
     h := g(u, s):  
     if verbose then print(u, h, `and`(op(is~(Im~(h) =~ 0.))) ); end if;
     if `and`(op(is~(Im~(h) =~ 0.))) then
        c := c, [h[], u]
     end if:
   end do:
   [c]
end proc:

# realCxyz(-1, 1, 1, 3, true)

PLOT3D( seq( CURVES(realCxyz(-2, 2, 0.005, s, false), COLOR(RGB, s/6, 0, 1-s/6), THICKNESS(s)), s=1..6), AXESLABELS('x(z)', 'y(z)', 'z'))

 

 


 

Download Example.mw

In a lot of languages the function to draw a graph names "plot" or "draw".
In Maple it's "plot", so I advice you to search the name "plot" in the help pages to get all the possibilitues of the function.

I don't understand what "relative" and "absolute" maxima mean to you.
So, look to the attached file, maybe it will help you.

 

Download Plot.mw

Here are Approximate and Exact poltynomial interpolations.

As I mentioned to tomleslie, tout data do not correspond to the plot your initial question contains


 

restart:

# OP's code after a few cosmetic arrangements

t := proc (x) options operator, arrow; b*x^10+c*x^9+d*x^8+e*x^7+f*x^6+g*x^5+h*x^4+i*x^3+j*x^2+k*x+l end proc;

solve([
        t(-2)  = 2,
        t(-1)  = 0,
        t(0)   = -2,
        t(1)   = 0,
        t(2.5) = -3,
        t(2.8) = 0,
        t(3.5) = 0,
        eval(diff(t(x), x), x = -2)  = 0,
        eval(diff(t(x), x), x = 0)   = 0,  
        eval(diff(t(x), x), x = 1)   = 0,
        eval(diff(t(x), x), x = 2.5) = 0
      ],
      [b, c, d, e, f, g, h, i, j, k, l]
);

plot(-.2688965311*x^10+.1687428810*x^9-.3166922031*x^8-1.490599337*x^7+1.885667707*x^6+3.845330330*x^5-6.939129440*x^4-2.523473874*x^3+7.112020606*x^2-2, x = -2..3.5):  # change of the range

proc (x) options operator, arrow; b*x^10+c*x^9+d*x^8+e*x^7+f*x^6+g*x^5+h*x^4+i*x^3+j*x^2+k*x+l end proc

 

[[b = -0.2688965311e-1, c = .1687428810, d = -0.3166922031e-1, e = -1.490599337, f = 1.885667707, g = 3.845330330, h = -6.939129440, i = -2.523473874, j = 7.112020606, k = 0., l = -2.]]

(1)

# My code

t  := proc (x) options operator, arrow; b*x^10+c*x^9+d*x^8+e*x^7+f*x^6+g*x^5+h*x^4+i*x^3+j*x^2+k*x+l end proc;
dt := unapply(diff(t(x), x), x);

proc (x) options operator, arrow; b*x^10+c*x^9+d*x^8+e*x^7+f*x^6+g*x^5+h*x^4+i*x^3+j*x^2+k*x+l end proc

 

proc (x) options operator, arrow; 10*b*x^9+9*c*x^8+8*d*x^7+7*e*x^6+6*f*x^5+5*g*x^4+4*h*x^3+3*i*x^2+2*j*x+k end proc

(2)

# data

# 1/ points

X  := [-2, -1,  0, 1, 2.5, 2.8, 3.5];
Y  := [ 2,  0, -2, 0,  -3,   0,   0];

# 2/ derivatives

dX := [-2, 0, 1, 2.5];
dY := [ 0, 0, 0,   0];

[-2, -1, 0, 1, 2.5, 2.8, 3.5]

 

[2, 0, -2, 0, -3, 0, 0]

 

[-2, 0, 1, 2.5]

 

[0, 0, 0, 0]

(3)

# sys  : system to solve
# vars : unknowns

sys  := [ (t~(X) =~ Y)[], (dt~(dX) =~ dY)[] ];
vars := convert(indets(t(1)), list);

[1024*b-512*c+256*d-128*e+64*f-32*g+16*h-8*i+4*j-2*k+l = 2, b-c+d-e+f-g+h-i+j-k+l = 0, l = -2, b+c+d+e+f+g+h+i+j+k+l = 0, 9536.743164*b+3814.697266*c+1525.878906*d+610.3515625*e+244.140625*f+97.65625*g+39.0625*h+15.625*i+6.25*j+2.5*k+l = -3, 29619.67667*b+10578.45595*c+3778.019983*d+1349.292851*e+481.890304*f+172.10368*g+61.4656*h+21.952*i+7.84*j+2.8*k+l = 0, 275854.7354*b+78815.63867*c+22518.75391*d+6433.929688*e+1838.265625*f+525.21875*g+150.0625*h+42.875*i+12.25*j+3.5*k+l = 0, -5120*b+2304*c-1024*d+448*e-192*f+80*g-32*h+12*i-4*j+k = 0, k = 0, 10*b+9*c+8*d+7*e+6*f+5*g+4*h+3*i+2*j+k = 0, 38146.97266*b+13732.91015*c+4882.812500*d+1708.984375*e+585.93750*f+195.3125*g+62.500*h+18.75*i+5.0*j+k = 0]

 

[b, c, d, e, f, g, h, i, j, k, l]

(4)

sol := solve(sys, vars)[] ;  # identical to the OP's

[b = -0.2688965311e-1, c = .1687428810, d = -0.3166922031e-1, e = -1.490599337, f = 1.885667707, g = 3.845330330, h = -6.939129440, i = -2.523473874, j = 7.112020606, k = 0., l = -2.]

(5)

# checking

T  := unapply(eval( t(x), sol), x):
dT := unapply(eval(dt(x), sol), x):

 T~( X) -~  Y;  # should be close to 0
dT~(dX) -~ dY;  # should be close to 0

[-0.9e-7, 0., 0., 0., -0.9e-7, 0.77e-6, 0.312e-5]

 

[0.19e-6, 0., 0., -0.2e-7]

(6)

plot(eval(t(x), sol), x=min(X[], dX[])..max(X[], dX[]))

 

# 1/ points

X  := [-2, -1,  0, 1, 5/2, 14/5, 7/2];
Y  := [ 2,  0, -2, 0,  -3,   0,   0];

# 2/ derivatives

dX := [-2, 0, 1, 5/2];
dY := [ 0, 0, 0,   0];

# 3/ exact solution

sys := [ (t~(X) =~ Y)[], (dt~(dX) =~ dY)[] ]:

sol := solve(sys, vars)[];

# 4/ checking

T  := unapply(eval( t(x), sol), x):
dT := unapply(eval(dt(x), sol), x):

 T~( X) -~  Y;
dT~(dX) -~ dY;

[-2, -1, 0, 1, 5/2, 14/5, 7/2]

 

[2, 0, -2, 0, -3, 0, 0]

 

[-2, 0, 1, 5/2]

 

[0, 0, 0, 0]

 

[b = -8786501831/326761419600, c = 551386670273/3267614196000, d = -206965577023/6535228392000, e = -19482815598053/13070456784000, f = 203690408213/108020304000, g = 50260227787849/13070456784000, h = -8245235930497/1188223344000, i = -4122869858861/1633807098000, j = 663981147407/93360405600, k = 0, l = -2]

 

[0, 0, 0, 0, 0, 0, 0]

 

[0, 0, 0, 0]

(7)

 


 

Download PolynomialFit.mw

If I understant correctly your problem, you want to assess the value of y2[n] = F(x2[n]) for each n=1..maxx.

If it is so here is a solution.
I have no doubt that someone else will prvide you a more concise/elegant solution


Sorry, I still get this ?!*&?!!  error when I try to load the content of the worksheet 
Maple Worksheet - Error
Failed to load the worksheet /maplenet/convert/Spline.mw .

 

Download Spline.mw

The simplest form of a third order equation is Z^3+p*Z+q = 0
This form can always be obtained after elementary transformations of the equation is a*Z^3+b*Z^2+c*Z+d = 0.

It's well known that the narure of the roots of the equation  Z^3+p*Z+q = 0 (either 3 pairwise distinct real roots, or 1 double real root plus a distinct simple one, or 1 real root plus 2 complex conjugate roots) depends on the sign of its discriminent R = 4*p^3+27*q^2.
Maple does not presents the solution this way, but it always present the three roots in complex form.
This means that these forms, which are general enough to cover the three situations R < 0, R = 0 and R > 0 are complicated, unless trivial situation (for instance q=0, or p=0 for instance).


 

solve(x^3+p*x+q, x)

(1/6)*(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)-2*p/(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3), -(1/12)*(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)+p/(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)+((1/2)*I)*3^(1/2)*((1/6)*(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)+2*p/(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)), -(1/12)*(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)+p/(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)-((1/2)*I)*3^(1/2)*((1/6)*(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)+2*p/(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3))

(1)

 


 

Download ThirdOrder.mw

Now do back to your problem.
One can show (see the attached file below), that the solutions with respect to B can be written B1=f(A) and B2=-f(A).
Injecting these solutions in a remaining equation "in A" results in a third order equation of the form a*A^3+b*A^2+A*x+d = 0.
It's then easy to identify the coefficients p and q of the standard form Z^3+p*Z+q = 0 (file below).

This means the "Z solutions" Z1, Z2, Z3 are of the form given in the file ThirdOrder.mw.
The "A solutions" are An,= Zn + b/(3*a).

The standard third order equation writes 
When solved for Z, the expression one obtains for Z1, Z2, Z3 are still rather simple.
So are the corresponding values of A1, A2, A3.

The difficulties come when the roots B1 and B2 are computed from f(An) (last line of the file below)

Unless you can specify yourself the sign of the discriminant R and do some work by hand, it doesn't seem possible to obtain simpler expressions of the solution (see also Edgardo's answer)

SimplifiedForm.mw

What is your need:

  • You want to run your code by clicking on some icon which would run a graphic user interface.
    Then use Maplets for instance
     
  • You want to run a .mw file without opening a Maple sesion before.
    Then use the online cmaple (windows) or maple (linux/OSX) (with all it's limitations concerning plots)
     
  • You want to run a .mw filr  from another one, just if this latter was a procedure.
    Use DocumentTools:-RunWorksheet
     
  • You want your code to be executed while opened in a worksheet.
    Look to the entry autoexecutedisabled ine the help pages
     
  •  any other...?

 

@Harry Garst 

Sorry for this late answer but I was busy  with some other stuff and I just forgot this thread.
Concerning the EM algorithm to identify the (known number of) components of a mixture of non necessarily gaussian distributions, I wrote something for personal use times ago.
This is just an implementation of the EM algorithm (classical one without jumps to account for situations with unknown number of components).
If you're interested I can deliver it after some modifications (and comments). But I'm not sure it will fit your educational purposes requirements.

The R code for tetra/poly choric correlations is rather long and complex to translate into Maple (package psych, function tetrachoric). I will ry to do something simpler, at least fro tetrachoric correlations.
I'll keep you in touch in the few days to come.

For now on here is a "raw" version of the EM algorithm I use

 

restart:

with(Statistics):

local D;
Digits := 20:
K  := 3:

 


STEP 1: Simulate a TrueMixture of 3 random variables

NbOfComponents := 3:

TrueMixture := table([
  # 1 = table([law = ExponentialDistribution, weight = 0.5, param = [4]       ]),
  # 1 = table([law = NonCentralChiSquare    , weight = 0.4, param = [3, 2]    ]),
  # 1 = table([law = MoyalDistribution      , weight = 0.4, param = [-1, 1]   ]),
    1 = table([law = NormalDistribution     , weight = 0.4, param = [0, 1]    ]),
    2 = table([law = NormalDistribution     , weight = 0.4, param = [3, 1]    ]),
    3 = table([law = NormalDistribution     , weight = 0.2, param = [10, 3]   ])
]);


for k from 1 to K do
   G__||k := RandomVariable(TrueMixture[k][law](op(TrueMixture[k][param]))):
end do:


# The sampling can be done in wide variety of ways
# The simpler, non purely stochastic method is draw a sample
# for each G__k proportional to its weight:
N    := 10^4:

P    := [seq(TrueMixture[k][weight], k=1..NbOfComponents)];

NS   := floor~(N *~ P);
S__1 := Sample(G__1, NS[1]):
S__2 := Sample(G__2, NS[2]):
S__3 := Sample(G__3, NS[3]):
S    := < seq((S__||k)^+, k=1..K) >;

Histogram(S):

for k from 1 to K do
   o__||k := TrueMixture[k][weight]:
   m__||k := op(1, TrueMixture[k][param]):
   s__||k := op(2, TrueMixture[k][param]):
end do:



 

TrueMixture := table([1 = table([law = NormalDistribution, weight = .4, param = [0, 1]]), 2 = table([law = NormalDistribution, weight = .4, param = [3, 1]]), 3 = table([law = NormalDistribution, weight = .2, param = [10, 3]])])

 

P := [.4, .4, .2]

 

NS := [4000, 4000, 2000]

 

Matrix(%id = 18446744078239458838)

(1)


STEP 2: Estimate the parameters of the mixture

mroll := rand(-5.0 .. 15.0);
sroll := rand(0 .. 10.0)

proc () options operator, arrow; RandomTools:-Generate(float('range' = -5.0 .. 15.0, 'method' = 'uniform')) end proc

 

proc () options operator, arrow; RandomTools:-Generate(float('range' = 0 .. 10.0, 'method' = 'uniform')) end proc

(2)

# Initialize Mixture with the correct distributions

MaxMix := 3:

Mixture := table([
  seq(k = table([law = TrueMixture[k][law], weight = 1/MaxMix, param = [mroll(), sroll()] ]), k=1..MaxMix)
]):



# EM algorithm

choix := combinat:-randperm([$(1..N)])[1..1000]:

verbose := false:


#-------------------------------------------------------------
# Initialisation
#
# (could be done differently)
for k from 1 to MaxMix do
   o__||k := Mixture[k][weight]:
   m__||k := Mixture[k][param][1]:
   s__||k := Mixture[k][param][2]:
end do:

if verbose then print("m = ", seq(m__||k, k=1..MaxMix)) end if;
if verbose then print("s = ", seq(s__||k, k=1..MaxMix)) end if;

 

#-------------------------------------------------------------
# algo.


printf(cat(seq("-------------------------", k=1..MaxMix), "---\n"));
printf(cat(seq(cat("   |     Component ", k, "     "), k=1..MaxMix), "  |\n"));
printf(cat(seq("-------------------------", k=1..MaxMix), "---\n"));
printf(cat(seq("   | weight  par 1  par2 ", k=1..MaxMix), "  |\n"));
printf(cat(seq("-------------------------", k=1..MaxMix), "---\n"));

R := 20:
for r from 1 to R do
   for k from 1 to MaxMix do
      pdf    := t -> Mixture[k][weight] * PDF(Mixture[k][law](op(Mixture[k][param])), t);
      L__||k := pdf~(S[choix]);
   end do:

   LL := `+`~(seq(L__||k, k=1..MaxMix));
   if verbose then print(seq(L__||k, k=1..MaxMix), LL) end if;
   
   for k from 1 to MaxMix do
      post__||k := L__||k /~ LL;
   end do:
   if verbose then print(seq(post__||k, k=1..MaxMix)) end if;

   for k from 1 to MaxMix-1 do
      o__||k := Mean(post__||k)[1]
   end do:
   o__||MaxMix := 1-add(o__||k, k=1..MaxMix-1):
   if verbose then print("o = ", seq(o__||k, k=1..MaxMix)) end if;

   for k from 1 to MaxMix do
      m__||k := add(post__||k *~ S[choix]) /~ add(post__||k);
   end do:
   if verbose then print("m = ", seq(m__||k, k=1..MaxMix)) end if;

   for k from 1 to MaxMix do
      s__||k := sqrt(add(post__||k *~ (S[choix] -~ m__||k)^~2) /~ add(post__||k));
   end do:
   if verbose then print("s = ", seq(s__||k, k=1..MaxMix)) end if;
   
   for k from 1 to MaxMix do
      if numelems(Mixture[k][param])=1 then
         Mixture[k][param] := m__||k;
      else
         Mixture[k][param] := [m__||k, s__||k];
      end if
   end do:

   if verbose then print();print() end if;
 
   printf(
           cat("%2d ", seq("  %-23a", k=1..MaxMix), "\n"),
           r,
           op(evalf[4]([seq([o__||k, op(Mixture[k][param])], k=1..MaxMix)]))
   )
end do:

------------------------------------------------------------------------------

   |     Component 1        |     Component 2        |     Component 3       |
------------------------------------------------------------------------------
   | weight  par 1  par2    | weight  par 1  par2    | weight  par 1  par2   |
------------------------------------------------------------------------------

 1   [.2836, 4.398, 4.942]    [.2630, 3.411, .8269]    [.4534, 1.707, 3.358]  

 2   [.3541, 4.668, 4.999]    [.2497, 3.125, .8237]    [.3962, 1.225, 2.419]  

 3   [.3373, 5.380, 4.960]    [.2579, 2.997, .8450]    [.4048, .8177, 1.760]  

 4   [.3169, 5.947, 4.793]    [.2726, 2.974, .8691]    [.4104, .5431, 1.444]  

 5   [.3028, 6.301, 4.652]    [.2919, 2.981, .8756]    [.4053, .3468, 1.261]  

 6   [.2936, 6.525, 4.552]    [.3078, 2.979, .8683]    [.3986, .2152, 1.141]  

 7   [.2885, 6.660, 4.482]    [.3192, 2.965, .8593]    [.3923, .1296, 1.064]  

 8   [.2862, 6.734, 4.436]    [.3269, 2.945, .8551]    [.3869, .7433e-1, 1.017]

 9   [.2853, 6.773, 4.406]    [.3322, 2.925, .8557]    [.3825, .3820e-1, .9876]

10   [.2850, 6.793, 4.387]    [.3360, 2.907, .8593]    [.3789, .1386e-1, .9693]

11   [.2851, 6.804, 4.375]    [.3389, 2.893, .8642]    [.3761, -.3123e-2, .9573]

12   [.2851, 6.810, 4.366]    [.3410, 2.881, .8692]    [.3738, -.1536e-1, .9492]

13   [.2852, 6.814, 4.360]    [.3427, 2.872, .8740]    [.3721, -.2442e-1, .9435]

14   [.2853, 6.817, 4.356]    [.3440, 2.864, .8781]    [.3707, -.3126e-1, .9394]

15   [.2853, 6.819, 4.353]    [.3451, 2.858, .8816]    [.3696, -.3649e-1, .9363]

16   [.2854, 6.821, 4.350]    [.3459, 2.854, .8846]    [.3687, -.4055e-1, .9340]

17   [.2854, 6.822, 4.348]    [.3466, 2.850, .8870]    [.3680, -.4371e-1, .9323]

18   [.2854, 6.823, 4.347]    [.3471, 2.847, .8889]    [.3675, -.4618e-1, .9309]

19   [.2854, 6.824, 4.345]    [.3476, 2.845, .8905]    [.3670, -.4813e-1, .9299]

20   [.2854, 6.825, 4.344]    [.3479, 2.843, .8918]    [.3667, -.4967e-1, .9290]

 

FormalMixture := evalf[4](add(o__||k * Mixture[k][law](op(evalf[4](Mixture[k][param]))), k=1..MaxMix));

.2854*NormalDistribution(6.825, 4.344)+.3479*NormalDistribution(2.843, .8918)+.3667*NormalDistribution(-0.4967e-1, .9290)

(3)

mix := add(Mixture[k][weight] * PDF(Mixture[k][law](op(Mixture[k][param])), t), k=1..MaxMix):
plots:-display(
   Histogram(S[choix]),
   plot(mix, t=min(S[choix])..max(S[choix]), color=red, thickness=3),
   seq(plot(Mixture[k][weight] * PDF(Mixture[k][law](op(Mixture[k][param])), t), t=min(S[choix])..max(S[choix]), color=gold), k=1..MaxMix)
)

 

 


 

Download EM_algorithm.mw

The error message indicates you cannot write 1/F__i = 0 .. 10  as a range.
Write instead F__i=0.1..SomeBigNumber.

Doing this you will get the clear error message
Error, (in Plot:-Inequality) additional variables not matching range variables found: {A__c, Z__b, Z__t, k__b, k__t, `(e__0)__mp`, M__max__B, M__min__B, sigma__ci, sigma__cs, sigma__ti, sigma__ts}
saying there exist non instanciated quantities in the inequalities to plot: just give them values

As suggested in the code below, use Explore to help you giving values tio the remaining indets
 

restart:

with(plots):

eqns :=
  {
    F__i*(1-e__0/k__b)/A__c+M__min__B/Z__t >= sigma__ti,
    .8*F__i*(1-e__0/k__t)/A__c-M__max__B/Z__b >= sigma__ts,
    .8*F__i*(1-e__0/k__b)/A__c+M__max__B/Z__t <= sigma__cs,
    F__i*(1-e__0/k__t)/A__c-M__min__B/Z__b <= sigma__ci, abs(e__0) <= `(e__0)__mp`
  };

inequal(eqns, F__i = 0.1..100, e__0 = -5 .. 5, color = "Niagara 2")

{sigma__ti <= F__i*(1-e__0/k__b)/A__c+M__min__B/Z__t, sigma__ts <= .8*F__i*(1-e__0/k__t)/A__c-M__max__B/Z__b, .8*F__i*(1-e__0/k__b)/A__c+M__max__B/Z__t <= sigma__cs, F__i*(1-e__0/k__t)/A__c-M__min__B/Z__b <= sigma__ci, abs(e__0) <= `(e__0)__mp`}

 

Error, (in Plot:-Inequality) additional variables not matching range variables found: {A__c, Z__b, Z__t, k__b, k__t, `(e__0)__mp`, M__max__B, M__min__B, sigma__ci, sigma__cs, sigma__ti, sigma__ts}

 

indets(eqns, name) minus {F__i, e__0}

{A__c, Z__b, Z__t, k__b, k__t, `(e__0)__mp`, M__max__B, M__min__B, sigma__ci, sigma__cs, sigma__ti, sigma__ts}

(1)

Explore(inequal(eqns, F__i = 0.1..100, e__0 = -5 .. 5, color = "Niagara 2"))

 


 

Download inequal.mw

 

By change of coordinates

Change_coordinates.mw
 

restart;

g := (x, y) -> min((3+((x-y)^(2))/(10)+(abs(x+y))/(sqrt(2))),(-abs(x-y)+(7)/(sqrt(2))));

proc (x, y) options operator, arrow; min(3+(1/10)*(x-y)^2+abs(x+y)/sqrt(2), -abs(x-y)+7/sqrt(2)) end proc

(1)

# the terms (x-y) and (x+y) suggest changing frame [0, x, y]into frame [O, u, v] by a rotation of angle Pi/2

xy2uv := [u=(x+y)/sqrt(2), v=(x-y)/sqrt(2)];
uv2xy := op(solve(xy2uv, [x, y]));

[u = (1/2)*(x+y)*2^(1/2), v = (1/2)*(x-y)*2^(1/2)]

 

[x = (1/2)*2^(1/2)*(u+v), y = (1/2)*2^(1/2)*(u-v)]

(2)

G := unapply(simplify(subs(uv2xy, g(x, y))), [u, v]);

proc (u, v) options operator, arrow; min(-(1/2)*2^(1/2)*(2*abs(v)-7), 3+(1/5)*v^2+abs(u)) end proc

(3)

q := 0;
H := (u, v) -> Heaviside(G(u, v)-q);  # means H(u, v) is either 0 or 1
 

0

 

proc (u, v) options operator, arrow; Heaviside(G(u, v)-q) end proc

(4)

L := 10:
plots:-contourplot(H(u, v), u=-L..L, v=-L..L, contours=[0.999], filled=true, grid=[100, 100]);  # blue = 1, red = 0

 

# let

f := x -> exp(-x^2/2)/sqrt(2*Pi);

# then the integrand is

phi := (x, y) -> 'h(x, y)' * f(x) * f(y):

# apply the rotation

Phi := (u, v) -> simplify('H(u, v)' * subs(uv2xy, f(x)*f(y))):

Phi(u, v)

proc (x) options operator, arrow; exp(-(1/2)*x^2)/sqrt(2*Pi) end proc

 

(1/2)*H(u, v)*exp(-(1/2)*u^2-(1/2)*v^2)/Pi

(5)

# then p := int(int(h(x, y)*exp((-x^2-y^2)*(1/2))/(2*Pi), y = -infinity .. infinity), x = -infinity .. infinity);
# just becomes  p := int(int(Phi(u, v), u = -infinity .. infinity), v = -a..a);
# where a is the half width of the blue strip.
#
# More of this, H(u, v)=1 inside this plur strip, then

p := int(int(f(u)*f(v), u = -infinity .. infinity), v = -a..a);  # here "a" is still unknown

erf((1/2)*2^(1/2)*a)

(6)

H := unapply(subs(uv2xy, h(x, y)), [u, v])

proc (u, v) options operator, arrow; h((1/2)*2^(1/2)*(u+v), (1/2)*2^(1/2)*(u-v)) end proc

(7)

# Find the position of the horizontal lines which bound the blue strip

solutions := [solve(G(u, v)=0, v)];

strip := sort(select(is, solutions, real));

eval(p, a=strip[2]);

evalf(%);

[7/2, -7/2, I*(5*abs(u)+15)^(1/2), -I*(5*abs(u)+15)^(1/2)]

 

[-7/2, 7/2]

 

erf((7/4)*2^(1/2))

 

.9995347418

(8)

 


 

Download Change_coordinates.mw
 

I'm going to bed, here are some elements to answer points 1 to 4.
5 is obviously ExcelTools:-Export(....)

I have no idea of your knowledge about all this stuff, so pardon me if reading this worksheet makes you feel  I'm condescending to you.

restart:

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

with(Statistics):

Simple model set the framework of Ordinary Least Squares (OLS)
X1 and X2 are deterministic variables
B1 and B2 are parameters of unknown values

E is a random Error (with dome adhoc properties [homoskedasticity]  I assume verified here [OLS framework again]
Y = B1*X1 +B2*X2 + E

Given a collection if triples (x1[n], x2[n], y[n]), find approximations b1 and b2 of B1 and B2


Step 1: simulation

N := 100:
SampleOfX__1 := Sample(Uniform(0, 10), N): # just to "create" x1 points in an arbitrary way
SampleOfX__2 := Sample(Uniform(0, 10), N): # just to "create" x1 points in an arbitrary way

Error := Sample(Normal(0, 10), N); # usual choice (= unbiased observation error + homoskedasticity)


B1 := 3 ; B2 := -4:   # remember these values are UNKNOWN in practical situations

ExactResponse := B1 *~ SampleOfX__1 +~ B2 *~ SampleOfX__2:

ObservedResponse := ExactResponse +~ Error:
 

Error := Vector(4, {(1) = ` 1 .. 100 `*Vector[row], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

3

(2)

# Rotate the plot to verify the points do not lie on a single plane (because of the additive error)

ScatterPlot3D(< SampleOfX__1, SampleOfX__2, ObservedResponse>^+, symbol=solidcircle, symbolsize=20);

 

Estimation of B1 and B2

OLS_result := LinearFit(
                         [x__1, x__2],
                         < SampleOfX__1, SampleOfX__2>^+,
                         ObservedResponse^+,
                         [x__1, x__2],
                         output=[parametervalues, leastsquaresfunction, residualstandarddeviation]
                       )

OLS_result := [Vector(2, {(1) = 3.24372606226451, (2) = -3.69397058902010}), 3.24372606226451*`#msub(mi("x"),mi("1"))`-3.69397058902010*`#msub(mi("x"),mi("2"))`, 10.2959448666143]

(3)

Interpretation:
The best linear unbiased estimation (BLUE) of B1 equals 3.24372...
The best linear unbiased estimation (BLUE) of B2 equal -3.39397...
The residual standard deviation has value 10.29...  
This values characterizes the "spreading" of the points around the optimal OLS plane base on the BLUE estimators above


 

plots:-display(
   ScatterPlot3D(< SampleOfX__1, SampleOfX__2, ObservedResponse>^+, symbol=solidcircle, symbolsize=20),
   plot3d(OLS_result[2], x__1=0..10, x__2=0..10, color=blue, transparency=0)  # adjust transparency and/or rotate
);

 

The main problem with the residualstandarddeviation is that it only gives an information about the quality of representation of the "BLUE" plane
Then this error is named "REPRESENTATION Error" or also "LEARNING Error"

What we are often interested in is to assess the representability of the model OUTSIDE of the points it was constructed on.
TRhat is to assess the GENERALIZATION error
The classical strategy to realize this is to use a resampling procedure.
For instance :
    1/ Leave One Out (LOO) where all the points but one are used to build the optimal model while the discarded one is used to assess the generalization error
        Repeat this for all the points and assemble these N estimators ot obtain the desired generalization error


   2/ Leave k Out (LkO) similar to LOO but with k points discarded (N!/k!/(N-k)! local estimators can obtained that way)

  3/  Bootstrap (simpler to refer to the corresponding Maple's help page)

  4/ Cross Validation, or CV, (close to LkO): split the DATA base (the N triples here) in a LEARNING base used to build the optimal model, plus a TEST or GENERALISATION
      base to assess its generalization error.
      As always repeat the process.

About the 25%-75% splitting parts I have a few things to say. Some others will use 1/3-2/3.
There is no strict rule beyond the fact that the model must remain  constructible on the choosen fraction of the DATA base
 

# Simulation

Data_base := < SampleOfX__1, SampleOfX__2, ObservedResponse>^+:

Learning_ratio := 0.25:
Learning_size  := round(Learning_ratio*N);
Test_size      := N - Learning_size;

KeptForLearning := combinat[randperm]([$1..N])[1..Learning_size];  # randomly select 25% of the data base
KeptForTest     := convert({$1..N} minus {KeptForLearning[]}, list);
 

25

 

75

 

[13, 1, 53, 86, 89, 59, 34, 75, 69, 42, 61, 90, 46, 74, 11, 44, 72, 35, 6, 84, 87, 27, 22, 45, 91]

 

[2, 3, 4, 5, 7, 8, 9, 10, 12, 14, 15, 16, 17, 18, 19, 20, 21, 23, 24, 25, 26, 28, 29, 30, 31, 32, 33, 36, 37, 38, 39, 40, 41, 43, 47, 48, 49, 50, 51, 52, 54, 55, 56, 57, 58, 60, 62, 63, 64, 65, 66, 67, 68, 70, 71, 73, 76, 77, 78, 79, 80, 81, 82, 83, 85, 88, 92, 93, 94, 95, 96, 97, 98, 99, 100]

(4)

Learning_base := Data_base[KeptForLearning, ..];
Test_base     := Data_base[KeptForTest, ..]    ;
 

Learning_base := Vector(4, {(1) = ` 25 x 3 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

Matrix(%id = 18446744078279553862)

(5)

# Learnt (=Trained) model

CV_result := LinearFit(
                         [x__1, x__2],
                         Learning_base,
                         [x__1, x__2],
                         output=[parametervalues, leastsquaresfunction, residualstandarddeviation]
                       )

CV_result := [Vector(2, {(1) = 3.52749864138740, (2) = -3.52762830240327}), 3.52749864138740*`#msub(mi("x"),mi("1"))`-3.52762830240327*`#msub(mi("x"),mi("2"))`, 9.05975305600152]

(6)

# Its predictions

CV_model := unapply(CV_result[2], [x__1, x__2]);

CV_predictions := Vector[column](Test_size, n -> op(CV_model(entries(Test_base[n, 1..2]))))

CV_model := proc (`#msub(mi("x"),mi("1"))`, `#msub(mi("x"),mi("2"))`) options operator, arrow; 3.5274986413873965*`#msub(mi("x"),mi("1"))`+(-1)*3.527628302403268*`#msub(mi("x"),mi("2"))` end proc

 

Vector[column](%id = 18446744078279555902)

(7)

# Its generalization error:

CV_Gen_Error := StandardDeviation(CV_predictions -~ Test_base[..,3])

HFloat(10.58269840059723)

(8)

A value not very far from the one of the representation error obtained in the OLS test, but it is related to the fact that the TRUE model is ALSO the
regression model (you could change the TRUE model by writting for instance Y = B1*X1 +B2*X2^2 + E and keep the regretion model Y = B1*X1 +B2*X2
unchanged: you will then see large differences between the trainin and the test errors)

Let's generate now R realizations of this test error:

R := 100: # it's just illustrative

CV_All := NULL:

for r from 1 to R do
  KeptForLearning := combinat[randperm]([$1..N])[1..Learning_size];  
  KeptForTest     := convert({$1..N} minus {KeptForLearning[]}, list);
  Learning_base   := Data_base[KeptForLearning, ..];
  Test_base       := Data_base[KeptForTest, ..]    ;
  CV_result       := LinearFit(
                         [x__1, x__2],
                         Learning_base,
                         [x__1, x__2],
                         output=[parametervalues, leastsquaresfunction, residualstandarddeviation]
                    ):
  CV_model        := unapply(CV_result[2], [x__1, x__2]);
  CV_predictions  := Vector[column](Test_size, n -> op(CV_model(entries(Test_base[n, 1..2])))):
  CV_Gen_Error    := StandardDeviation(CV_predictions -~ Test_base[..,3]):
  
  CV_All := CV_All, [CV_result[], CV_Gen_Error]:
end do:

CV_All := [CV_All]:

All_Gen_Errors := op~(-1, CV_All);
 

[HFloat(9.651751409842975), HFloat(11.073297744209476), HFloat(10.582065315202824), HFloat(11.370905480893054), HFloat(10.66062548594883), HFloat(10.197860333787814), HFloat(12.41984696724287), HFloat(10.26697182510627), HFloat(11.282596469419106), HFloat(10.953890045733575), HFloat(10.805611200949137), HFloat(10.269181816780955), HFloat(11.738688054559878), HFloat(11.038255433740698), HFloat(10.288690840528693), HFloat(10.302616068603378), HFloat(9.92059362228101), HFloat(11.369276018689709), HFloat(10.531187095376493), HFloat(10.528157324620775), HFloat(10.794603635544576), HFloat(10.687373377037863), HFloat(10.093956416633349), HFloat(9.765740103319256), HFloat(10.04186222004179), HFloat(10.960507825341455), HFloat(10.681157271733277), HFloat(10.95331271884746), HFloat(10.49331163470941), HFloat(10.606985758324798), HFloat(10.322822452553421), HFloat(10.831019288718446), HFloat(10.690102166263637), HFloat(10.707222219250296), HFloat(11.496121117213429), HFloat(10.202813494517128), HFloat(10.870415933896236), HFloat(10.228559285899076), HFloat(9.901316980202198), HFloat(11.235324777668138), HFloat(10.902603625122918), HFloat(10.914117298801584), HFloat(11.456526480997686), HFloat(10.09331224647256), HFloat(11.537103244138848), HFloat(10.13310322832185), HFloat(10.928632140543822), HFloat(10.206451530904193), HFloat(11.120120136270378), HFloat(9.729935966292539), HFloat(10.186651165055293), HFloat(10.098925551423877), HFloat(11.25433646374303), HFloat(10.736132209751101), HFloat(11.02680833986016), HFloat(10.109821507353962), HFloat(10.426438708344806), HFloat(10.283779712420099), HFloat(10.695264381144266), HFloat(10.31222432561016), HFloat(10.89859930478856), HFloat(11.080714909982953), HFloat(10.153367981091353), HFloat(10.322124849198017), HFloat(9.635392933723857), HFloat(9.67503021569015), HFloat(10.990485781584994), HFloat(10.580001065747831), HFloat(10.98300811348783), HFloat(10.819589530631744), HFloat(10.530170012359214), HFloat(10.980387844133084), HFloat(11.454124102451425), HFloat(11.917027702421562), HFloat(10.560728962174448), HFloat(10.303471444788272), HFloat(10.247638749777725), HFloat(10.53435615420061), HFloat(10.74631482502725), HFloat(11.091629780383638), HFloat(10.435095538197611), HFloat(10.592305219979297), HFloat(10.128061684406104), HFloat(10.343475968239149), HFloat(10.593868706318588), HFloat(9.677543534261027), HFloat(10.375985506158676), HFloat(9.944500315401031), HFloat(10.774385240991624), HFloat(10.274926237145262), HFloat(10.967918055475566), HFloat(10.74535182383358), HFloat(10.861251445658647), HFloat(10.471201895181126), HFloat(9.6799192328998), HFloat(11.025356009868046), HFloat(10.519396354811088), HFloat(10.4657467263297), HFloat(11.225574737612337), HFloat(10.16744656966274)]

(9)

# the better estimate of the generalization error of a bilinear model trained on 25% of the data base is

GENERALIZED_Error := sqrt(Mean(All_Gen_Errors^~2));

HFloat(10.630183338660515)

(10)

# Matrix of all the BLUE b1 and b2
#
# Scatter plot of all the BLUEs from all the test bases, plus (red asterixk) the OLS BLUE

CV_BLUE := convert(op~(1, CV_All), Matrix)^+;

plots:-display(
   plot(CV_BLUE, style=point, color=black),
   plot([[entries(OLS_result[1], nolist)]], color=red, symbol=asterisk, symbolsize=30, style=point)
);

CV_BLUE := Vector(4, {(1) = ` 100 x 2 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

 

# Let's built a very simple aggregated model

Aggregated_Model := Mean(CV_BLUE[..,1]) * x__1 + Mean(CV_BLUE[..,2]) * x__2;

HFloat(3.295250552411322)*x__1-HFloat(3.7572438843065274)*x__2

(11)

In some sense, this aggregated model can be said to have a generalization error equal to   
Typesetting:-mrow(Typesetting:-mi("GENERALIZED_Error", italic = "true", mathvariant = "italic"), Typesetting:-mo("&coloneq;", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.2777778em", rspace = "0.2777778em"), Typesetting:-mn("10.6301833386605", mathvariant = "normal"), Typesetting:-mspace(height = "0.0ex", width = "0.0em", depth = "0.0ex", linebreak = "newline"), Typesetting:-mi("From", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("this", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo("error", bold = "false", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("point", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("of", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("view", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("it", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("is", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("a", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("better", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("model", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("than", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("the", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("original", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("OLS", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("model", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo("for", bold = "false", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("which", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("NO", bold = "true", italic = "false", mathvariant = "bold", fontweight = "bold"), Typesetting:-mo(" ", bold = "true", mathvariant = "bold", fontweight = "bold", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("generalization", bold = "true", italic = "false", mathvariant = "bold", fontweight = "bold"), Typesetting:-mo(" ", bold = "true", mathvariant = "bold", fontweight = "bold", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo("error", bold = "true", mathvariant = "bold", fontweight = "bold", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo(" ", bold = "true", mathvariant = "bold", fontweight = "bold", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("can", bold = "true", italic = "false", mathvariant = "bold", fontweight = "bold"), Typesetting:-mo(" ", bold = "true", mathvariant = "bold", fontweight = "bold", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("be", bold = "true", italic = "false", mathvariant = "bold", fontweight = "bold"), Typesetting:-mo(" ", bold = "true", mathvariant = "bold", fontweight = "bold", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("assessed", bold = "true", italic = "false", mathvariant = "bold", fontweight = "bold"), Typesetting:-mspace(height = "0.0ex", width = "0.0em", depth = "0.0ex", linebreak = "newline"), Typesetting:-mo(" ", bold = "true", mathvariant = "bold", fontweight = "bold", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mspace(height = "0.0ex", width = "0.0em", depth = "0.0ex", linebreak = "newline"), Typesetting:-mo(" ", bold = "true", mathvariant = "bold", fontweight = "bold", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("But", italic = "false", mathvariant = "normal"), Typesetting:-mo(",", mathvariant = "normal", fence = "false", separator = "true", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.3333333em"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo("for", bold = "true", mathvariant = "bold", fontweight = "bold", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("the", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("representation", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("point", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("of", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("wiew", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("this", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("aggregated", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("model", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("is", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("worse", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("than", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("the", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("OLS", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("model", italic = "false", mathvariant = "normal"), Typesetting:-mo(":", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.2777778em", rspace = "0.2777778em"), Typesetting:-mspace(height = "0.0ex", width = "0.0em", depth = "0.0ex", linebreak = "newline"), Typesetting:-mi("The", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mi(""), Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi("N", italic = "true", mathvariant = "italic"), Typesetting:-mo("&minus;", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.2222222em", rspace = "0.2222222em"), Typesetting:-mn("1", mathvariant = "normal")), mathvariant = "normal")), Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi("N", italic = "true", mathvariant = "italic"), Typesetting:-mo("&minus;", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.2222222em", rspace = "0.2222222em"), Typesetting:-mn("3", mathvariant = "normal")), mathvariant = "normal")), linethickness = "1", denomalign = "center", numalign = "center", bevelled = "false"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("factor", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("comes", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo("from", bold = "true", mathvariant = "bold", fontweight = "bold", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("technical", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("considerations", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("ense", italic = "false", mathvariant = "normal"), Typesetting:-mo(",", mathvariant = "normal", fence = "false", separator = "true", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.3333333em"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi("N", italic = "false", mathvariant = "normal"), Typesetting:-mo("&minus;", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.2222222em", rspace = "0.2222222em"), Typesetting:-mn("3", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("degrees", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("of", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("freedom", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo("to", bold = "false", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("assess", italic = "false", mathvariant = "normal"), Typesetting:-mo(" ", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.0em", rspace = "0.0em"), Typesetting:-mi("residualstandarddeviation", italic = "false", mathvariant = "normal")), mathvariant = "normal"), Typesetting:-mspace(height = "0.0ex", width = "0.0em", depth = "0.0ex", linebreak = "newline"), Typesetting:-mspace(height = "0.0ex", width = "0.0em", depth = "0.0ex", linebreak = "auto"))

Predictons_of_Aggregated_Model := Vector[column](N, n -> op(CV_model(entries(Data_base[n, 1..2]))));

Representation_error_of_Aggregated_Model := StandardDeviation(Predictons_of_Aggregated_Model -~ Data_base[.., 3]) * (N-1) / (N-3)

Predictons_of_Aggregated_Model := Vector(4, {(1) = ` 1 .. 100 `*Vector[column], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

HFloat(10.490314893535418)

(12)

 


 

Download OLS.mw

 

Kitonum has been  faster than me !!!

I was prepared to delete my reply, very close to Kitonum's with whom I agree on a lot of points....
excepted the final result which, for a reason I do not know about, is not the same .
So I finally decided to send "my" proposal... maybe there is an obvious typo error in it or there is something Kitonum maybe will be able to explain?


 

restart

k := T -> 2/cosh(2/T)/coth(2/T)

proc (T) options operator, arrow; 2/(cosh(2/T)*coth(2/T)) end proc

(1)

local delta, e:

delta := T -> sqrt(1-k(T)^2*sin(theta)^2)

proc (T) options operator, arrow; sqrt(1-k(T)^2*sin(theta)^2) end proc

(2)

e := T -> -2*tanh(2/T)*diff(k(T)/2/Pi*Int(sin(theta)^2/delta(T)/(1+delta(T)), theta=0..Pi), T)

proc (T) options operator, arrow; -2*tanh(2/T)*(diff((1/2)*k(T)*(Int(sin(theta)^2/(delta(T)*(1+delta(T))), theta = 0 .. Pi))/Pi, T)) end proc

(3)

ee := evalf(e(T)):  # e(T) does the differentiation accortinng to T and evalf replaces Int by int

eval(ee, T=1):      # substitute T by 1 in ee

evalf(%);           # compute the result

-.3586637209

(4)

 


 

Download proposal.mw

You made a confusion between 
c[1] := 2; c[2] := 5; 
and farther
c := (beta[1]*s(t)*p(t).(lambda[2](t)-lambda[1](t)))/(w[2]*(a+p(t)))+(lambda[2](t).e(t)+(i(t)+alpha.e(t)).lambda[3](t)-(gamma.i(t)+p(t)).lambda[4](t))/w[2]; 

Then, in the lambda[2](t) edo, the first term of the rhs contains c[1] and appears as (beta[1]*s(t)*p(t)......./w[2])instead of (I guess) 2.
Idem in the lambda[2](t) edo, the first term of the rhs contains c[2] and appears as (beta[1]*s(t)*p(t)......./w[2])instead of (I guess) 5.
This explains the error you get Error, (in fproc) unable to store '-1.*HFloat(0.0)[1]': the HFloat(0.0)[1] refers to (beta[1]*s(t)*p(t)......./w[2])


It's to you to correct this.
For example by writting 
C:= (beta[1]*s(t)*p(t).(lambda[2](t)-lambda[1](t)))/(w[2]*(a+p(t)))+(lambda[2](t).e(t)+(i(t)+alpha.e(t)).lambda[3](t)-(gamma.i(t)+p(t)).lambda[4](t))/w[2]; 
u[2] := min(max(0, C), 1);


Once done it should be work better, but don't be in a hurry, with maxmesh=2400 Maple takes its time to provide a solution (about 5 minutes).

Download MaybeThis.mw

Is this suit you?
(don't worry about the expression of A: I just copied-pasted the content of your question).

PS: Given a Graph G, the command AdjacencyMatrix(G) returns the equivalent of your matrix A, but, I do not know about a Maple's procedure which construct a graph G given its AdjacencyMatrix.
 

restart

with(GraphTheory):

A := Matrix(4, 10, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 1, (2, 5) = 0, (2, 6) = 1, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 1, (3, 1) = 1, (3, 2) = 1, (3, 3) = 0, (3, 4) = 0, (3, 5) = 1, (3, 6) = 0, (3, 7) = 1, (3, 8) = 1, (3, 9) = 1, (3, 10) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 0})

A := Matrix(4, 10, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 1, (2, 5) = 0, (2, 6) = 1, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 1, (3, 1) = 1, (3, 2) = 1, (3, 3) = 0, (3, 4) = 0, (3, 5) = 1, (3, 6) = 0, (3, 7) = 1, (3, 8) = 1, (3, 9) = 1, (3, 10) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 0})

(1)

edges := map(x -> if A[op(x)] = 1 then {x[]} end if, {indices(A)});

{{1, 3}, {2, 3}, {2, 4}, {2, 6}, {2, 10}, {3, 5}, {3, 7}, {3, 8}, {3, 9}}

(2)

G := Graph(edges);

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], Array(%id = 18446744078096242318), `GRAPHLN/table/1`, 0)

(3)

DrawGraph(G, style=spring);

 

 


 

Download Graph.mw

@a_simsim 

A third way is "table" ; maybe easier to manipulate (for instance to add or to remove a "field")


Table.mw

PS: a "field" can itself be a table, and so on

Personally I prefer using records

I find them more readable than arrays (Tom's proposal) and the structure is less artificial than the one obtained with arrays.
Betond, I don't have enough enough skill in Maple to tell you which of these two structures is more efficient in terms of memory usage or manipulation time.

Last by not least, as you mention Matlab, record is closer than struct than array is
 

restart:

#              the
#              name
#              you
#              want
#               |
#               |
#               V
r := Record(
             'MyArray' = Array(1..3, 1..3,1..3),
             'data'    = Vector([2,1,3,5,6]),
             'status'  = true,
             'x'       = 123.456,
             'n'       = 17,
             'name'    = "astring",
             'MyVect'  = Vector([true, false, true, true, true, false]),
             'MyMat'   = Matrix(7,7, symbol=x),
             'MyPlot'  = plot(x^2, x=0..5)
            ):

r:-data

Vector[column]([[2], [1], [3], [5], [6]])

(1)

r:-n

17

(2)

r:-MyPlot

 

 


 

Download Record.mw

1 2 3 4 5 6 7 Page 3 of 9