@acer 

A way (the only one?) to place the text in the foreground and the hatching in the background.

Text_foreground_hatches_background.mw

@janhardo 

Download Improvements.mw

@Carl Love

I don't think I wrote something like that.

As a reminder I wrote

Then @Carl Love probably confused by your relation (4) told you you were computed a sum and Maple  equation (5) was right.

and

So I'm sorry to say that if you really want to compute mixtures properties, all the answers that have been given here contain, are wrong in a general sense, largely due to the confusion you initiated between the weighted sum and the mixture  of random  variables.

OP's question began with him writting
I would like to build a Gaussian Mixture Model in Maple using the most straightforward approach:
and providing what I took as an illustration.
The original error comes from the OP defining a mixture of two gaussian random variables as a sum of two gaussians random variables, which is generally false unless the two (more generally all) random variables have the sam expectation) which led you in a direction that does not correspond to the one the use of the word mixture should imply.

You write
I'd like to point out the I didn't say a single word about mixtures of random variables, nor did I claim to
indeed, for the reason I explained above you focused on OP's equation (4) not on the mixture word.

I think you misinterpreted my response and overreacted to my words, as it was not my intention to offend you.
I consider that the misunderstanding is over.

@janhardo 

The formula for the pdf you write in mixtureDist is correct, so is the one for the mean, but your variance formula is correct only if the two gaussian random variables share the same mean.

For decades, I have used C++ and FORTRAN debuggers.
They have continuously improved over the years, enabling debugging parallel codes on massively parallel architectures while becoming more and more friendly and easy to use.
Maple's debugger is really pitiful in comparison.

The only situation I keep using it is to access built-in procedures: when I have to debug my own code I prefer doing it the old-fashioned way, that is using print/printf commands instead.

@CarlosRodrigues

I thought that a week would have been enough time for you to analyze the answers you received and have the courtesy to respond to us.

@acer 

as indeed a good idea.
Even if this took me a lot of time to arrive where I wanted to, O finally discovered that Statistics:-LinearFit and Statistics:-Fit both use at the very end the NAG routine g02daf.
So debugging endeds here but I got nevertheless some interesting informations.

Thanks for your help

@acer 

I don't know what to say.

Indeed this proc calls LinearFitting:-LinearLS (which I can't find any trace of), likely some local names, when I expected a code detailing the way the fitted model is computed.

In a few words, I just noticed accidentally that Statistics:-LinearFit (and Statistics:-Fit) didn't give the expected result on a quite simple test case, was issuing a meaningless warning, and that the help pages contained incorrect information.

When I write "expected result" I mean the result one can obtain "by hand" using different solving strategies, results which are the same than the one LinearAlgebra:-LeastSquares  or (an hammer to kill a fly) Optimization:-LSSolve both provide.
After digging into the subject it appears that Statistics:-LinearFit probably uses its own code, that the solution method it truly uses is not the one the help page says it uses (which "Normally" [I quote the help page] is the QR decomposition), that it likely imposes its own Digits value (for the warning iti displayed can be avoided increasing Digits), ... and a lot of other things.

With you reply it looks like the snake which bites its own tail and I'm a bit in limbo

I was writting a post on the subject and expected a definitive answer to conclude this post and send it.
Maybe filling a RCS too where I would demand that  Statistics:-LinearFit uses LinearAlgebra:-LeastSquares (whose help page contain reliable informations).

If you want (and if you have time to spend on this subject) I can provide you with a worksheet which details all of this.

@nm 

The first thing is not to compute EllipticE(x, m) but EllipticE(sin(x), m)  (note this will work only where the application x → sin(x) is a one-to-one map given the way EllipticE is defined in help(EllipticE) ).
The reason why is given in the attached file, see equation (6).

Doing this goves MMA_EllipticE(x, m) = Maple_Elliptic(sin(x), sqrt(m) ) iif -Pi/2 <= x <= Pi/2.

To get the same results for any x value you must to "periodize" EllipticE this way:

• For instance, to get the same result than MMA_EllipticE(3, -1), you must compute 
  Maple_Elliptic(sin(Pi/2), I ) + Maple_Elliptic(3-sin(Pi/2), I ).
  The procedure EllipticE_a_la_MMA shows you how to do this.
   
• ... or you use function MA in equation (2) below, which seems too me far simpler.

Download EllipticE__MMA_vs_Maple.mw

@nm  @dharr

Maple notation and MMA/SageMath notation differ.

help(EllipticK)

Elliptic integrals and the related functions are well described in the Table of Integrals Series and Products, Gradshteyn and Ryzhik (G&R) and in the popular Handbook of Mathematical Functions edited by Abramowitz and Stegun (A&S). In A&S, these functions are expressed in terms of a parameter m, representing the square of the modulus k entering the definition of the Elliptic, JacobiPQ and InverseJacobiPQ functions in Maple and G&R. For example, the  K(m) function shown in A&S is numerically equal to the Maple EllipticK(sqrt(m)) command.

Both  MMA and SageMath use the Abramowitz-Stegun notation.
From §17.3 p 590
So Maple returning EllipticE(sinx(x), I) is the same thing than MMA or SageMath retyrning EllipticE(sin(x), -1)

@janhardo 

Of course Maple 2015 should have ALWAYS return an error when computing R := A.B.A^+

The problem is it didn't  when A and B were floats.

@C_R 

I did a mistake... but Maple 2015 is gugged in a way which balances my mistake.
Look at this: A being a vector column Maple 2015 computes cvorrectly A.B.A⁺ if A and B hve floating point elements!

Download Maple_2015_Bug.mw

In my procedure I mistakenly wrote Wres  := (res . Rinv . res^+); where res is a vector column.
With the Maple 2015's bug above I get the correct result.
In your more recent version (I guess) this bug has been fixed, hence the error you are receiving..

Here is my corrected worksheet Iterative_Gaussian_Process_Emulator.mw
Apologies

@C_R 

In your comments you use L as a synomym for lambda. Correct? : right, sorry for my not having been consistent in my notations.

About your error:

1. In Maple 2015 X0, Y0, X and Y are column vectors: is it so with your Maple Version?
   1. If not convert all of them into Vector[column],
   2. if yes then : 
      1. Download the attached file and duplicate it.
      2. Open them both in your Maple version.
      3. Keep the first one as the reference (do not execute it).
      4. Execute the second one.
      5. Tell me where the error occurs.

Compare_outputs_2015_your_version.mw

@C_R

In particular your informations about TuringBot are very interesting and I will havea look on it soon.

In case you would be interested in doing the same kind of thing in the future, here is a procedure which builds a GPE (Gaussian process Emulator) of an unknown scalar function f : X → Y (X being 1D) iteratively.
There is no stopping criterion: the decision to stop the process is yours as soon as you consider the GPE fits the data points well enough.

The process starts by selecting a small size subset of the data set (which may be interesting if this latter is large)  and by constructing a first GPE based on this subset.
Next the process enters an iterative step where a new point (belonging to the data set) is added to the subset at each iteration and the GPE is rebuilt.
The selection of the new point is based on the prediction error of the GPE built at the previous iteration.

The whole process converges quite rapidly.
For your data on can get a very good result by using less than 15% of the data set.

Iterative_Gaussian_Process_Emulator.mw

@dharr 

The attached file contains a variant named QR_RationalFractionFit of procedure RationalFractionFit .
Compared to this latter I replaced the call to LinearFit by a the QRDecomposition of the matrix named inputs plus a resolution of the system R . C = Q⁺ . output  using LinearSolve.

The attached file presents the [9, 9] case, which generates a Warning, model is not of full rank. You will see that the  QR decomposition stabilizes the system an avoid the occurrenceof this warning.

The new QR_RationalFractionFit procedure is incomplete because it does not return the many statistics that LinearfFit returns. So there is still a lot of work to be done butnothing complicated here, simply some time to spend on it.

QR_RationalFractionFit.mw




Using this  QR_RationalFractionFit procedure enables to compute the AIC criterion up to values greater than those at which you had to stop because ot the warnings.
The best model from the AIC criterion point of view has a numerator of degree dn = 11 and a denominator of degree dd =  13, corresponding to a 23-parameter model.
Comparing this number to the 101 points the data set contains means that this "best" model has a complexity of about 23% and so that it does not overfit the data.



Note that the fittings have been done up to the [20, 20] rational fraction without receiving any warning.