You will see that the circular curve disappears in the last frame. This has to do with rounding errors. The last sqrt is a complex number with a small imaginary part. You can avoid this by replacing the third command by

C := animatecurve( [1/2*p,1/2*Re(sqrt(4-p^2)), p = 0..2], frames=50 ):

@GrecoRoman Firstly, in Error.mw the formal parameter is data, not x. Secondly, the solve command in EuclideanRegression has multiple solutions as can be seen by adding a print command:

Moreover, in the last unapply statement the y must be replaced by x

In Euclidean.mw the third solution is apparantly the solution you want.

@GrecoRoman Firstly, in Error.mw the formal parameter is data, not x. Secondly, the solve command in EuclideanRegression has multiple solutions as can be seen by adding a print command:

Moreover, in the last unapply statement the y must be replaced by x

In Euclidean.mw the third solution is apparantly the solution you want.

@kjetil87  , this is the right approach.
The sum command tries to calculate a symbolic sum if the range contains parameters or infinity. If the range contains only numbers you can better use the add command. When calculating a symbolic sum, only the generic values of the variables are taken into account; eventual singularities you have to find out for yourself. In your example it is easily seen that only p1=0 and p1=0 lead to special cases that have to do with "0⁰ = 1" So the best thing you can do is indeed

restart;
A := p1^abs(l)*p2^abs(l-2)*exp(-I*2*Pi*f*l*T):
sum( A, l=1..infinity); # generic case
                               p1             
                  ----------------------------
                  p2 (exp(2 I Pi f T) - p1 p2)
S := add(A, l=1..2) + sum(A, l=3..infinity): # special cases
eval(S,p2=0);
                        2                 
                      p1  exp(-4 I Pi f T)
eval(S,p1=0);
                               0

@kjetil87  , this is the right approach.
The sum command tries to calculate a symbolic sum if the range contains parameters or infinity. If the range contains only numbers you can better use the add command. When calculating a symbolic sum, only the generic values of the variables are taken into account; eventual singularities you have to find out for yourself. In your example it is easily seen that only p1=0 and p1=0 lead to special cases that have to do with "0⁰ = 1" So the best thing you can do is indeed

restart;
A := p1^abs(l)*p2^abs(l-2)*exp(-I*2*Pi*f*l*T):
sum( A, l=1..infinity); # generic case
                               p1             
                  ----------------------------
                  p2 (exp(2 I Pi f T) - p1 p2)
S := add(A, l=1..2) + sum(A, l=3..infinity): # special cases
eval(S,p2=0);
                        2                 
                      p1  exp(-4 I Pi f T)
eval(S,p1=0);
                               0

@jam to test if the left hand side of an equation equals the right hand side you can use ?evalb (after simplify).

@jam to test if the left hand side of an equation equals the right hand side you can use ?evalb (after simplify).

@Christopher2222 yes (see also Doug Meade's suggestion and Kitonum's answer).
The function that you want to plot is only defined on the circular disk x² + y² ≤ 25.

Therefore you can use polar coordinates for x and y to obtain constant ranges for the parameters:

plot3d( subs( {x=r*cos(t),y=r*sin(t)}, [x,y,sqrt(25-x^2-y^2)] ), r=0..5, t=0..2*Pi, axes=boxed );

Another possibility is

 plot3d( sqrt(25-x^2-y^2), x=-5..5, y=-sqrt(25-x^2)..sqrt(25-x^2), axes=boxed );

But this also misses a few points.

@Christopher2222 yes (see also Doug Meade's suggestion and Kitonum's answer).
The function that you want to plot is only defined on the circular disk x² + y² ≤ 25.

Therefore you can use polar coordinates for x and y to obtain constant ranges for the parameters:

plot3d( subs( {x=r*cos(t),y=r*sin(t)}, [x,y,sqrt(25-x^2-y^2)] ), r=0..5, t=0..2*Pi, axes=boxed );

Another possibility is

 plot3d( sqrt(25-x^2-y^2), x=-5..5, y=-sqrt(25-x^2)..sqrt(25-x^2), axes=boxed );

But this also misses a few points.

@acer , I also tried to find an indexing function, but I didn't succeed.
But your quick solution doesn't work:

restart;
m := LinearAlgebra:-RandomMatrix(4,6):
M := ArrayTools:-UpperTriangle(m);
V := Vector( convert( select( x -> x<>0, M ), list ) ):
G:=(m,n,W) -> Matrix(m,n,(i,j)->`if`(i>j,0,W[(i-1)*n-((i-1)*(i)/2)+j])):
G(4,6,V);

@acer , I also tried to find an indexing function, but I didn't succeed.
But your quick solution doesn't work:

restart;
m := LinearAlgebra:-RandomMatrix(4,6):
M := ArrayTools:-UpperTriangle(m);
V := Vector( convert( select( x -> x<>0, M ), list ) ):
G:=(m,n,W) -> Matrix(m,n,(i,j)->`if`(i>j,0,W[(i-1)*n-((i-1)*(i)/2)+j])):
G(4,6,V);

Casper, I'm not sure if I have all the indices right. But I think this may help you further.

restart:
para := proc(C,new)
   global p;
    # new = 0 or 1
    seq(log(p[j,C]/(1-p[j,C]))=mu+tau[j]+eta[C] + new*tau[j]*eta[C],j=2..4)
end proc:

getpVec:=proc(k,c,new)
  local i,j, parC;
  global p, wp;
  for i to c do
    solve({para(i,new)},[seq(p[j,i],j=2..(k+1))]);
    parC[i] := Vector(k, [rhs~(op(%))] );
  end do;
  add(wp[i]*parC[i],i=1..(c-1))+(1-add(wp[i],i=1..(c-1)))*parC[c];
end proc:
getpVec(3,3,1);

Casper, I'm not sure if I have all the indices right. But I think this may help you further.

restart:
para := proc(C,new)
   global p;
    # new = 0 or 1
    seq(log(p[j,C]/(1-p[j,C]))=mu+tau[j]+eta[C] + new*tau[j]*eta[C],j=2..4)
end proc:

getpVec:=proc(k,c,new)
  local i,j, parC;
  global p, wp;
  for i to c do
    solve({para(i,new)},[seq(p[j,i],j=2..(k+1))]);
    parC[i] := Vector(k, [rhs~(op(%))] );
  end do;
  add(wp[i]*parC[i],i=1..(c-1))+(1-add(wp[i],i=1..(c-1)))*parC[c];
end proc:
getpVec(3,3,1);

...so I make a procedure to do this:

RestoreMatrix := proc(r,c,V)
local M, i,j,k;
M := Matrix(r,c):
k := 1:
for i to r do
   for j to i do
     M[j,i] := V[k]:
     k := k+1
   end do
end do;
for i from r+1 to c do
   for j to r do
     M[j,i] := V[k];
     k := k+1
   end do
end do;
return M
end proc;

RestoreMatrix(4,6,V);

...so I make a procedure to do this:

RestoreMatrix := proc(r,c,V)
local M, i,j,k;
M := Matrix(r,c):
k := 1:
for i to r do
   for j to i do
     M[j,i] := V[k]:
     k := k+1
   end do
end do;
for i from r+1 to c do
   for j to r do
     M[j,i] := V[k];
     k := k+1
   end do
end do;
return M
end proc;

RestoreMatrix(4,6,V);