6 years, 74 days

@LichengZhang   @Carl LoveCarl...

@LichengZhang   @Carl Love

Carl Love already gave an elegant answer to the question "How many cycles of length n does a given Graph contain?" (which obiously includes your first one " I'd like to know  whether it contains a C4 (cycle of 4) as its subgraph"

I authorise myself to reproduce Carl's code here

 > restart:
 > with(GraphTheory): with(RandomGraphs):
 > randomize():
 > NV := 8: NE := 12: G  := RandomGraph(NV, NE):
 > # The code below has been written by Carl Love # https://www.mapleprimes.com/questions/226523-How-To-Determine-The-Number-Of-Circles-Of--Graph-#answer257002 kCycles:= (A::Matrix, n::posint, k::And(posint, Not({1,2})))->    add(        (-1)^(k-i)*binomial(n-i,n-k) *          add(LinearAlgebra:-Trace(A[S,S]^k), S= combinat:-choose(n,i)),        i= 2..k    )/2/k:
 > # How to use it (cf above link) A:= GraphTheory:-AdjacencyMatrix(G): seq(printf("G contains %2d cycles of length %2d\n", kCycles(A, NV, k), k), k= 3..NV);
 G contains  2 cycles of length  3 G contains  5 cycles of length  4 G contains  6 cycles of length  5 G contains  5 cycles of length  6 G contains  4 cycles of length  7 G contains  0 cycles of length  8
 >

Download Does_G_contains_a_given_cycle.mw

Does a given Graph contain K(n) where K(n) is the complete graph with n vertices?
GraphTheory:-CliqueNumber(G) returns the highest value of N fpr wich K(N) is a complete subgraph of G.
Then, G contains at least one complete subgraph K(n) for n=3..N.
Conversely, G contains no complete subgraph K(n) such that n > N

Here is an example (adapt it with  ...

Here is an example (your own file is treated later)

Be carefull: as recently mentioned in another thread, Maple doesn't support loops (FromVertex = ToVertex)

PS: if your graph is not oriented, you can replace
edges := { convert(M, listlist)[] }
by
edges := { convert~(convert(M, listlist), set)[] }

 > restart:
 > kernelopts(homedir): file := cat(%, "/Desktop/graph.xlsx"):
 > M := Import(file);
 (1)
 > M := map(round, M);
 (2)
 > edges := { convert(M, listlist)[] }
 (3)
 > G :=GraphTheory:-Graph(edges)
 (4)
 > GraphTheory:-DrawGraph(G)
 >

Download Excel2Graph.mw

Now the worksheet with your own file (2896 Vertices graph)
You will have a lot of work to do to obtain a readable display !

Excel2Graph_bis.mw

Here are plots for 1 and 2 categorical v...

Here are plots for 1 and 2 categorical variables.
They were done with Maple 2015. I do not have Maple 2018 on this laptop and I do not remember if Maple 2018 has, or not,  specific  features to plot categorical variables

For an unknown reason I can't load the content of the worksheet

CategoricalVariables.mw

PS to customize the 3D plots you might try other commands than
plottools:-line( [i, j, 0], [i, j, Sizes[i, j]], thickness=10, color=GenderColor[i]),
for instance, having set GenderColor := [blue, red]:
plottools:-cylinder([i, j, 0], 0.2, Sizes[i, j], color=GenderColor[i], style=surface),

Why not solving  symbolically ...

Why not solving  symbolically with formal parameters instead of numeric ones, and next substituting the formal parameters by the suitable numerical values?

 > ode1 := (2/3)*diff(theta(t), t\$2)+A__1*theta(t)-A__2*x(t) = 0;
 (1)
 > ode2 := 2*diff(x(t),t\$2)+B__1*x(t)-B__2*theta(t) =B__3*sin(B__4*t);
 (2)
 > ics := theta(0) = 0, (D(theta))(0) = D__1, x(0) = 0, (D(x))(0) = 0
 (3)
 > solution := dsolve([ode1, ode2, ics]):
 > simplify~(eval(solution, [A__1=1250, A__2=500, B__1=1000, B__2=500, B__3=50, B__4=20, D__1=1/100]))
 (4)
 > simplify~(eval(solution, [A__1=1250, A__2=500, B__1=1000, B__2=500, B__3=50, B__4=90, D__1=1/100]))
 (5)
 >

Download TryThis.mw

Here an answer obviously less concise th...

Here an answer obviously less concise than vv's answer, but I found funny to adress the problem from a pure geometrical point of view.

I had already done this exercise with Geogebra... and I wondered if Maple could do the same ?

Still this f... error
Maple Worksheet - Error
Failed to load the worksheet /maplenet/convert/geometry.mw .

Download geometry.mw

The problem is such general, with suchma...

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);
 (1)
 > sys := [ seq(randpoly(vars, degree = 2)=0, n=1..3) ]: print~(sys):
 (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:
 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]
 (3)
 > # Here are numerical values of the couple (x, y)  if z=0 (for instance) evalf~(subs~(z=0, rhs~(ExactSolutions[1][1..2])))
 (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(%);
 (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 dr...

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

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
 (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);
 (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];
 (3)
 > # sys  : system to solve # vars : unknowns sys  := [ (t~(X) =~ Y)[], (dt~(dX) =~ dY)[] ]; vars := convert(indets(t(1)), list);
 (4)
 > sol := solve(sys, vars)[] ;  # identical to the OP's
 (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
 (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;
 (7)
 >

Download PolynomialFit.mw

If I understant correctly your problem, ...

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

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

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

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

@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:
 (1)

STEP 2: Estimate the parameters of the mixture

 > mroll := rand(-5.0 .. 15.0); sroll := rand(0 .. 10.0)
 (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));
 (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 w...

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")
 > indets(eqns, name) minus {F__i, e__0}
 (1)
 > Explore(inequal(eqns, F__i = 0.1..100, e__0 = -5 .. 5, color = "Niagara 2"))
 >

Download inequal.mw

By change of coordinatesChange_coordinat...

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))));
 (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]));
 (2)
 > G := unapply(simplify(subs(uv2xy, g(x, y))), [u, v]);
 (3)
 > q := 0; H := (u, v) -> Heaviside(G(u, v)-q);  # means H(u, v) is either 0 or 1
 (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)
 (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
 (6)
 > H := unapply(subs(uv2xy, h(x, y)), [u, v])
 (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(%);
 (8)
 >

I'm going to bed, here are some elem...

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)
 (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:
 (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]                        )
 (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);
 (4)
 > Learning_base := Data_base[KeptForLearning, ..]; Test_base     := Data_base[KeptForTest, ..]    ;
 (5)
 > # Learnt (=Trained) model CV_result := LinearFit(                          [x__1, x__2],                          Learning_base,                          [x__1, x__2],                          output=[parametervalues, leastsquaresfunction, residualstandarddeviation]                        )
 (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]))))
 (7)
 > # Its generalization error: CV_Gen_Error := StandardDeviation(CV_predictions -~ Test_base[..,3])
 (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);
 (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));
 (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) );
 > # Let's built a very simple aggregated model Aggregated_Model := Mean(CV_BLUE[..,1]) * x__1 + Mean(CV_BLUE[..,2]) * x__2;
 (11)

In some sense, this aggregated model can be said to have a generalization error equal to

 > 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)
 (12)
 >

Download OLS.mw

 First 22 23 24 25 26 27 28 Page 24 of 30
﻿