mmcdara

You probably need to add assumptions...

Answered: mmcdara 7891

April 13 2023

0 1

Example

>

restart; with(LinearAlgebra)

>

TPE := (1/2)*(int(int(D__11*(diff(w, x, x))^2+2*D__12*(diff(w, x, x))*(diff(w, y, y))+4*D__66*(diff(w, x, y))^2+D__22*(diff(w, y, y))^2-2*q__0*w, x = 0 .. a), y = 0 .. b));

Warning, computation interrupted

>

J := D__11*(diff(w, x, x))^2+2*D__12*(diff(w, x, x))*(diff(w, y, y))+4*D__66*(diff(w, x, y))^2+D__22*(diff(w, y, y))^2-2*q__0*w:

>	J := simplify(J, size):

>	int(J, x=0..a, y=0..b) assuming a > 0, b > 0, i::posint

(1)

>

Download total_PE_with_assumptions.mw

Don't know if it may help ...

Answered: mmcdara 7891

April 13 2023

3 6

Before several customization features have been introduced in the GraphTheory package, I used to do my own customization.
It can be lengthy but quite simple if you observe the structure that DrawGraph displays.

I don't know what Maple 2023 is capable to do, but in case it can do exactly what you want, there is an example of a method to change the shape, position and style of the arrows or the style of the vertices.
procedures MoveArrow and ChangeVertex have to be used after you have exploited all the possibilities offered by GraphTheory.

>

restart:

>	interface(version)

(1)

>	with(GraphTheory): with(plots):

>	G := Digraph({[1, 2], [2, 3], [3, 4], [4, 1]}):

>	dg := DrawGraph(G):

>	print~([op(dg)]):

(2)

CHANGING EDGES

>

# The last POLYGONS structures describe each edges with its arrow.
#
# If you want to change the style and position of any arrow you can do this.

MoveArrow := proc(edge, pc, angle, length, closed, col)
  local f, dx, dy, alpha, beta1, beta2, head, end1, end2, arrow, style:
  f      := (dx, dy) -> piecewise(dx >= 0, arctan(dy/dx), Pi+arctan(dy/dx)):
  dx, dy := (op(1, edge)[2] -~ op(1, edge)[1])[];
  alpha := f(dx, dy);
  beta1 := alpha-angle;
  beta2 := alpha+angle;
  head   := op(1, edge)[1]*~pc +~ op(1, edge)[2]*~(1-pc):
  end1   := head -~ length *~ [cos, sin](beta1):
  end2   := head -~ length *~ [cos, sin](beta2):
  style := remove(type, [op(edge)], list)[]:
  if closed then
    arrow := [end1, head, end2, end1]:
   #return POLYGONS([op(1, edge)[], ListTools:-Reverse(op(1, edge))[]], style),
   #       POLYGONS(arrow, COLOR(RGB, op(ColorTools:-NameToRGB24(col))))
    return POLYGONS(arrow, COLOR(RGB, op(ColorTools:-NameToRGB24(col)))),
           POLYGONS([op(1, edge)[], ListTools:-Reverse(op(1, edge))[]], style)
  else
    arrow := [end1, head], [head, end2]:
    return POLYGONS(op(1, edge), arrow, style)
  end if:
end proc:

NV     := numelems(Vertices(G)):
NE     := numelems(Edges(G)):
pol    := select(has, [op(dg)], POLYGONS):
edges := pol[-NE..-1]:
others := pol[1], remove(has, [op(dg)], POLYGONS)[]:

cols     :=["Red", "Cyan","Yellow", "Magenta"]:
pos      := [0.4, 0.3, 0.4, 0.3]:
NewEdges := seq(MoveArrow(edges[n], pos[n], Pi/12, 0.1, true, cols[n]), n=1..NE):

display(PLOT(others, NewEdges), axes=none)

CHANGING VERTICES

>	centers := map(n -> Statistics:-Mean(convert(op(n, others[1]), Matrix)), [$1..NV]): diameters := map(n -> abs~(op(n, others[1])[1] -~ centers[n]), [$1..NV]):

diameters := [Vector[row](2, {(1) = HFloat(0.0288), (2) = HFloat(0.028799999999999937)}), Vector[row](2, {(1) = HFloat(0.0288), (2) = HFloat(0.0288)}), Vector[row](2, {(1) = HFloat(0.028799999999999937), (2) = HFloat(0.028799999999999937)}), Vector[row](2, {(1) = HFloat(0.028799999999999937), (2) = HFloat(0.0288)})]

(3)

>

ChangeVertex := proc(shape, c, d, theta, col)
  local fig, n:
  if shape = "disk" then
    fig := display(plottools:-disk([0, 0], d[1], color=col));
  elif shape = "ellipse" then
    fig := display(plottools:-ellipse([0, 0], entries(d, nolist), color=col, filled=true))
  elif shape[1] = "P" then
    n   := parse(shape[2..-1]):
    fig := display(polygonplot([seq([cos(2*Pi*i/n), sin(2*Pi*i/n)], i = 1..n)], color=col)):
    fig := plottools:-scale(fig, entries(d, nolist))
  end if;
  plottools:-translate(plottools:-rotate(fig, theta), entries(c, nolist));
end proc:

>

# shapes := {"disk", "ellipse", ...) (see plottools)
#           or "Pn" where n is an integer >= 3

NewDiameters := map(d -> d*~[2, 1], diameters):
rotations    := (Pi/NV)*~[$1..NV]:
shapes       := ["disk", "ellipse", "P4", "P7"]:
display(
  seq(ChangeVertex(shapes[n], centers[n], NewDiameters[n], rotations[n], "Chartreuse"), n=1..NV)
  , PLOT(others[2..-1], NewEdges)
  , axes=none
  , scaling=constrained
)

>

Download Customization_1.mw

I don't have time right now to answe...

Answered: mmcdara 7891

April 13 2023

1 5

I don't have time right now to answer your second question.
More of this it requires some precisions: for instance what are S[1] and S[2] for they do notappear in the expresison of Omega?

First question:
You have made an error in reasoning: see the explanation in the attached file after the yellow highlighted text.
Basically the variance and covariance are sufficient statistics from which other statistics can be build (the linear correlation coefficient for instance), not the opposite.
Thus you cannot do some transformation on random variables, believe the correlation coefficient is an invariant, and deduce the covariance from its value and the two variances.

RVs_sum_mmcdara.mw

About sufficient statistics see here https://en.wikipedia.org/wiki/Sufficient_statistic

Rescaling is not the very problem...

Answered: mmcdara 7891

April 11 2023

2 0

I agree with @dharr that it isalways better to rescale the data to avoid numerical problem.
But rescaling hides a major bug in Fit/LinearFit: both these functions say C1=0 while it's not.

I already observed several weaknesses or failures of LinearFit, this is just a new one.
In case you doubt the results LinearFit returns, use the formula to get the vector A of the regressor coefficients:

A := (Z^+ . Z)^(-1) . Z^+ . Y;

In the attached file you will see that the wrong result returned by Fit/LinearFit comes from a wrong estimation of the rnk of the information matrix ((Z^+ . Z) in the text): as this rank is rqual to 3 Fit/LinearFit finds it equal to 2 only for precision problems.

These issues can be solve by rescaling the data,; but rescaling them BECAUSE C1=0 is not a good argument.

>

restart:

>	with(Statistics):

>	interface(version)

(1)

>

Data := Matrix([[4.74593554708566, 11385427.62, 2735660038000], [4.58252830679671, 25469809.77, 12833885700000], [4.29311160501838, 1079325200, 11411813200000000], [4.24176959154225, 1428647556, 18918585950000000], [5.17263072694618, 1428647556, 18918585950000000], [4.39351114955735, 1877950416, 30746202150000000], [4.39599006758777, 1428647556, 18918585950000000], [5.79317412396815, 2448320309, 49065217290000000], [4.48293612651735, 2448320309, 49065217290000000], [4.19990181982522, 2448320309, 49065217290000000], [5.73518217699046, 1856333905, 30648714900000000], [4.67943831980476, 3071210420, 75995866910000000], [4.215240105336, 2390089264, 48670072110000000], [4.41566877563247, 3049877383, 75854074610000000], [4.77780395369828, 2910469403, 74061327950000000], [4.96617430604669, 1416936352, 18891734280000000], [4.36131111330988, 1416936352, 18891734280000000], [5.17783192063198, 1079325200, 11411813200000000], [4.998266287191, 1067513353, 11402362980000000], [4.23366152474871, 2389517120, 48661380410000000], [4.58252830679671, 758079709.3, 5636151969000000], [6.82390874094432, 1304393838, 14240754750000000], [4.24176959154225, 912963601.2, 8621914602000000], [4.52432881167557, 573965555.4, 3535351888000000], [4.84133601918601, 573965555.4, 3535351888000000], [6.88605664769316, 732571773.2, 5558875538000000], [4.35575841415627, 1203944381, 13430693320000000], [4.42527441640593, 955277678, 8795128298000000], [6.82390874094432, 997591754.9, 8968341995000000], [4.35144484433733, 143039477.1, 305355143300000]]):

>	X := Data[.., [2, 3]]: Y := Data[.., 1]: Fit(C1+C2v+C3w, X, Y, [v, w]); LinearFit(C1+C2v+C3w, X, Y, [v, w]);

Warning, model is not of full rank

(2)

>	Z := `<\|>`(Vector(numelems(Y), 1), X): r := proc(d) Digits:=d: printf("Digits %2d rank = %d\n", d, LinearAlgebra:-Rank(Z^+ . Z)); Digits:=10 end proc: r(10), r(20), r(40): Digits(); A := (Z^+ . Z)^(-1) . Z^+ . Y; F := `<,>`(1, x1, x2) . A assuming x1::real, x2::real;

Digits 10 rank = 1
Digits 20 rank = 2
Digits 40 rank = 3

A := Vector(3, {(1) = 4.659353816079307, (2) = 0.5985084089529947e-9, (3) = -0.27350964718426345e-16}, datatype = float[8])

(3)

>	plots:-display( ScatterPlot3D(`<\|>`(X, Y), symbol=solidsphere, color=blue, symbolsize=20), plot3d(F, x1=(min..max)(X[..,1]), x2=(min..max)(X[..,2]), style=wireframe, color=gray) );

WITH SCALED DATA

>	model := unapply( LinearFit(C1+C2v+C3w, Scale(X), Scale(Y), [v, w]), [v, w])

(4)

>	mX, sX := (Mean, StandardDeviation)(X); mY, sY := (Mean, StandardDeviation)(Y);

(5)

>	MODEL := model((x1-mX[1])/sX[1], (x2-mX[2])/sX[2]) * sY + mY

(6)

>	plots:-display( ScatterPlot3D(`<\|>`(X, Y), symbol=solidsphere, color=blue, symbolsize=20), plot3d(F, x1=(min..max)(X[..,1]), x2=(min..max)(X[..,2]), style=wireframe, color=gray), plot3d(MODEL, x1=(min..max)(X[..,1]), x2=(min..max)(X[..,2]), style=wireframe, color=red) );

>

Download Bug_In_Fit.mw

The explicit computation of A and LinearFit applied to Scaled data give the same result.
At the minimum the fact that LinearFit applied on raw data gives a wrong result is a weakness, but I lean towards a bug:

either LinearFit must automatically rescale the data,
eiither LinearFit must say "The information matrix is numerically not of full rank, so I don't return any result".

Note: the only possibility for this matrix not to be of full rank (3) is that there exist three numbers a, b, c such that.

a*X[..,1]+b*X[..,2]+c = 0

It's easy to verify this is not true:

ScatterPlot(X[..,1], X[..,2]);
# the residual sum of squares = 0 iif a*x1+b*x2+c = 0
LinearFit([1, x1], X[..,1], X[..,2], [x1], output=residualsumofsquares)
                     1.03756169900198*10^33
# other way
LinearAlgebra:-LinearSolve(Z, Vector(numelems(Y), 0))
               < a, b, c > =  < 0, 0, 0 >

It's simpler and less sensitive to roundoff errors to use this later command to check if the matrix (Z^+ . Z) is of full rank, or nor.
Note that you can use identically

LinearAlgebra:-LinearSolve(X, Vector(numelems(Y), 0))

At least compare the results these 4 commands return

LinearAlgebra:-Rank(Z^+ . Z);         # uncorrect result (should be 3)
LinearAlgebra:-Rank(Z);               # uncorrect result (should be 3)
                               1
                               2
LinearAlgebra:-Rank(X^+ . X);         # uncorrect result (should be 1)
LinearAlgebra:-Rank(X);               # correct result
                               1
                               2

This does not add much to what has been ...

Answered: mmcdara 7891

April 11 2023

3 0

This does not add much to what has been said before.
I just imagined that the challenge is to find "by hand" a formal development of P_n(x) ????

So the idea I followed is this one:

Start from n = 1 and interpret your equality this way: P₁(x) is the quotient of (x-x^₂)⁴ by (1+x²) and -4 is its remainder.
The division of (x-x^₂)⁴ by (1+x²) can easily be done by hand to find the expression of P₁(x) and the remainder -4.
The remainder of the division of (x-x^₂)⁴ by (1+x²) is obviously (-4)ⁿ .
The binomial expansion then gives the expression of P_n(x) as a function of P₁(x) and r₁=-4.

>

restart

>	f := n -> (x-x^2)^(4*n); g := (1+x^2);

(1)

>	# You have this relation f(n) = g*p[n] + (-4)^n

(2)

>	# An interpretation is "the ramainder of the division of f(n) by # (1+x^2) is (-4)^n and its quptient is (1-x)*p(n)(x): # # For n = 1: do the euclidian division of f(1) by (1+x^2) QR := proc(n) quo(f(n), (1+x^2), x, 'r'): return %, r end proc: p[1], r[1] := QR(1)

(3)

>	# then 'f(1)' = 'p[1]'*(1+x^2) + r[1]

(4)

>	# For n > 1 f(n) writestat 'f(n)' = ('p[1]'*(1+x^2) + 'r[1]')^n

(5)

>	# Thus the remainder of division of f(n) by (1+x^2) # is r[n] := r[1]^n

(6)

>	# And its quotient q[n] verifies ('p[1]'(1+x^2) + 'r[1]')^n = (1+x^2)p[n] + 'r[1]'^n

(7)

>	# Using the binomial expansion of the lhs from k=1 to k=n # (to remove the constant term) one gets the expression of # q[n] = P(n): P := n -> sum(factorial(n)/factorial(k)/factorial(n-k)(1+x^2)^(k-1)p[1]^k*r[1]^(n-k), k=1..n)

proc (n) options operator, arrow; sum(factorial(n)*(x^2+1)^(k-1)*p[1]^k*r[1]^(n-k)/(factorial(k)*factorial(n-k)), k = 1 .. n) end proc

(8)

>

P(1);

(9)

>	expand(P(2))

(10)

>

Download By_hand.mw

A way...

Answered: mmcdara 7891

April 09 2023

2 1

This is not exactly what you want, but it's not that far.
You can still tune this code to improce the rendering, in particulat the positions and fonts of the three parameter counters.

>

restart;

>	with(DEtools):with(plots):with(plottools):

>	sigma1:=e1alpha: sigma2:=e2delta:

>	g:=x/(1+betax^2); f:=(thetax-1)(1-x)(1+betax^2)-y; h:=alphag-z-sigma1; j:=delta*y-sigma2;

alpha*x/(beta*x^2+1)-e1*alpha-z

(1)

>

anim := proc()
  local p0, p1, p3, q0, q1, q2, q3:

  p0:=theta->display(
               plot([1/theta,y,y=0..1],linestyle=dash,color=green),
               textplot([0, 1.05, `#mo("θ")`=evalf[3](theta)], align=right)
             ):
  p1:=e1->display(
               plot([x,e1,x=0....1.5],color=blue),
               textplot([1.5, 1.05, `#mo("e1")`=evalf[3](e1)], align=left)
             ):
  p3:=beta->display(
               plot([x,x/(1+beta*x^2),x=0..1.5],color=magenta),
               textplot([0.75, 1.05, `#mo("β")`=evalf[3](beta)])
             ):

  q0:=animate(p0,[theta],theta=2...10):
  q1:=animate(p1,[e1],e1=0.1..1):
  q3:=animate(p3,[beta],beta=0..1.5):
  q2:=plot([1,y,y=0..1],linestyle=dash,color=green):

  [q0, q1, q2, q3]
end proc:

display(anim(),title=" ", view=[0..1.5,0..1.1])

>

Download animate_mmcdara.mw

For fun...maybe?...

Answered: mmcdara 7891

April 09 2023

0 4

Not cuneiform symbols but Tibetan marks.

Download For_fun.mw

Here is a list of some special ASCII characters: http://www.addressmunger.com/special_ascii_characters/
Maple can visualize (almost) all of these characters up to

𝐀

Capital A bold; Capital 𝐀 bold

𝐀

excluded

plots:-intersectplot...

Answered: mmcdara 7891

April 09 2023

2 0

But this requires to build the cartesian equation of the two tori:

with(plottools):
with(plots):

# torus([1,1,1], 1, 2)
# torus([1,6,1], 1, 2)

# Torus with symmetry plane // (O,x, yj) 
#   center       = [p, q, r], 
#   major radius = a, 
#   minor radius = b.
#
# 
#
# Cartesian equation

ce := (p, q, r, a, b) ->  ((x-p)^2 + (y-q)^2 + (z-r)^2 + a^2 - b^2)^2 = 4*a^2*((x-p)^2 + (y-q)^2):
  
# Cartesian equations of torus 1 and torus 2
T1 := ce(1, 1, 1, 2, 1);
T2 := ce(1, 6, 1, 2, 1);

# Their intersection is contained in the plane y=7/2

plots:-intersectplot(T1, T2, x=-2..4, y=3.4..3.6, z=0..2, color=red, thickness=5, grid=[20, 2, 20])

Download Intersection.mw

restart:with(Statistics):X := Vector([1,...

Answered: mmcdara 7891

April 09 2023

2 0

restart:
with(Statistics):
X := Vector([1, 2, 3, 4, 5, 6], datatype=float):
Y := Vector([2, 5.6, 8.2, 20.5, 40.0, 95.0], datatype=float):
f := unapply(PowerFit(X, Y, v), v)
                                 HFloat(2.041674366793877)
v -> HFloat(1.4795777879958958) v                         
f(Pi)
                   HFloat(15.316373717844828)

Why do you use a deprecated function?Not...

Answered: mmcdara 7891

April 08 2023

0 3

Why do you use a deprecated function?
Note that using solve instead doesn't change things, but at least you avoid unnecessary problems.

The fact is that the right and left limit of your expression at x=1/2 are not the same. More precisely they real parts are thesame, but they have opposite imaginary parts.

>

restart:

>	interface(version)

(1)

>	expr := (x-1/2)(1-sqrt(1+3x(x-2/3)/(5(x-1/2)^2)))/(x-2/3): sol := solve(expr);

(2)

>	limit(expr, x=sol[1], left); limit(expr, x=sol[1], right); convert(series(expr, x=1/2, 2), polynom);

((3/5)*I)*5^(1/2)*csgn(I/(x-1/2))+(-6+((12/5)*I)*5^(1/2)*csgn(I/(x-1/2)))*(x-1/2)

(3)

>	limit(expr, x=sol[2], left); limit(expr, x=sol[2], right);

(4)

>	plots:-dualaxisplot( plot(Re(expr),x=-1/4..3/4, color=blue, numpoints=400, thickness=3, discont=true, legend="Re"), plot(Im(expr),x=-1/4..3/4, color=red, numpoints=400, thickness=3, discont=true, legend="Im"), gridlines=true )

>

Download solve_Re_Im.mw

A lot of errors:...

Answered: mmcdara 7891

April 08 2023

0 3

Within the loop: remove the square brackets in the rhs of the two expressions.
Within the loop: replace sum by add.
Within the loop: you wrote (bracket removed)
```
F*(sum(lambda^i*u[i], i = 0 .. n))
```
instead of
```
F(sum(lambda^i*u[i], i = 0 .. n))
```
Note that F being undefined
You missed to define u[0] before entering the loop:
```
u[0] := 1+t*e^(x+y)
```
The definition of A[n], at least in Maple 2015.2, is not correct when n=0
The integral is not correctly written:
1. dt*dt means nothing
2. int(..., t=0..t, t=0..t) is not correct, (Maple 2015.2), write instead
```
int(int(..., t=0..t) t=0..t) 
```
... and I missed probably some other mistakes.

Here is a code which works, but I do not pretend it implements correctly the Adomian's method as you define it (browse this site, a lot of questions about it have already been posted and answered).
Try and analyze it and correctly while respecting Maple's syntax.

>

restart

>	PDEtools[declare](prime = x):

(1)

>	# The two next lines are useless # equ1 := u[tt] = 1/2(u[xx]+u[yy])+u^2: # ICS := u(x, y, 0) = 1, u[t](x, y, 0) = e^(x+y), lambda = 0, u[0] = 1+te^(x+y), F(u[0]) = u[0]^2

>

F    := v -> v^2:
u[0] := 1+t*e^(x+y);
A[0] := F(u[0]);

for n from 0 to 5 do
  u[n+1] := (1/2)*int(int(diff(add(u[i], i = 0 .. n), x, x)+diff(add(u[n], i = 0 .. n), y, y), t=0 .. t), t=0..t)
            +
            int(A[n], t = 0 .. t):
  A[n+1] := diff(F(sum(lambda^i*u[i], i = 0 .. n+1)), lambda$(n+1))/factorial(n);
end do:

(2)

>	# Extremely lengthy expressions

Download ADM-1_mmcdara.mw

Lasr point: even using simplify(..., size) the expressions become rapidly hugely complex.
So do not try do visualize u[n] for n > 3.

How to correlate N gaussian random varia...

Answered: mmcdara 7891

April 07 2023

2 7

As you have all the formulas for step2 and beyond, I answer only your first question.
The attached file shows how to build a random vector whose gaussian components are linearly correlated.

The underlined bold text defines the limits of the method: all random variables have to be gaussian and the correlation must ne linear.
So do not try to extend this to other situations, such as, for instance, the linear correlation of uniform random variables: the result will be false.

>

restart

>	with(LinearAlgebra): with(Statistics): with(RealDomain):

>	N := 3: MarginalStDev := Vector(N, symbol=sigma); cor_Matrix := Matrix(N$2, (i, j) -> `if`(j=i, 1, rho[i, j]), shape=symmetric); cov_Matrix := DiagonalMatrix(MarginalStDev) . cor_Matrix . DiagonalMatrix(MarginalStDev);

MarginalStDev := Vector(3, {(1) = sigma[1], (2) = sigma[2], (3) = sigma[3]})

$cor_Matrix := Matrix(3, 3, {(1, 1) = 1, (1, 2) = rho[1, 2], (1, 3) = rho[1, 3], (2, 2) = 1, (2, 3) = rho[2, 3], (3, 3) = 1}, storage = triangular[upper], shape = [symmetric])$

$cov_Matrix := Matrix(3, 3, {(1, 1) = sigma[1]^2, (1, 2) = sigma[1]*rho[1, 2]*sigma[2], (1, 3) = sigma[1]*rho[1, 3]*sigma[3], (2, 1) = sigma[1]*rho[1, 2]*sigma[2], (2, 2) = sigma[2]^2, (2, 3) = sigma[2]*rho[2, 3]*sigma[3], (3, 1) = sigma[1]*rho[1, 3]*sigma[3], (3, 2) = sigma[2]*rho[2, 3]*sigma[3], (3, 3) = sigma[3]^2})$

(1)

>

# Conditions for cor_Matrix to be positive definite

rhos := indets(cor_Matrix):

# Sylvester's criterion: the N main minors must be strictly positive.

det_11 := cor_Matrix[1, 1]: # this minor isobviously positive
det_22 := Minor(cor_Matrix, 3, 3):
det_33 := Determinant(cor_Matrix):

rho_conds := solve({det_22 > 0, det_33 > 0}, rhos):
print~(rho_conds):

(2)

>	# Cholesky decomposition assuming rho_conds are met conditions := map(x -> op([x > -1, x < 1]), indets(cor_Matrix))[], map(x -> x > 0, indets(MarginalStDev))[]; LUDecomposition(cov_Matrix, method='Cholesky'): L := simplify(%) assuming conditions;

Warning, Matrix is non-hermitian; using only data from upper triangle

Parse:-ConvertTo1D, "arguments to %1 must be (integer,posint)", "_Inert_RATIONAL"

(3)

>	# Correlation by linear combination: # # step 1: Let A a random variable of independent centered gaussian variables A := `<,>`(seq(RandomVariable(Normal(0, 1)), i=1..N));

$A := Vector(3, {(1) = _R, (2) = _R0, (3) = _R1})$

(4)

>	# step 2: Define the random vector B this way B := L . A

B := 3

(5)

>	# step 3: check B 'Expectation(B)' = Mean~(B); 'Cov(B)' = simplify(Mean~(B . Transpose(B)))

"?=RTABLE(18446744074197627966,MATRIX([[0], [0], [0]]),Vector[column])"

$Cov(B) = (Matrix(3, 3, {(1, 1) = sigma[1]^2, (1, 2) = sigma[1]*rho[1, 2]*sigma[2], (1, 3) = sigma[1]*rho[1, 3]*sigma[3], (2, 1) = sigma[1]*rho[1, 2]*sigma[2], (2, 2) = sigma[2]^2, (2, 3) = sigma[2]*rho[2, 3]*sigma[3], (3, 1) = sigma[1]*rho[1, 3]*sigma[3], (3, 2) = sigma[2]*rho[2, 3]*sigma[3], (3, 3) = sigma[3]^2}))$

(6)

>	# final step: Define the random vector C as B + the vector of desired means. # Check C C := B + Vector(N, symbol=mu): 'Expectation(C)' = Mean~(C); 'Cov(C)' = simplify( Mean~(C . Transpose(C)) - Mean~(C) . Transpose(Mean~(C)) )

"?=RTABLE(18446744074197628238,MATRIX([[mu[1]], [mu[2]], [mu[3]]]),Vector[column])"

$Cov(C) = (Matrix(3, 3, {(1, 1) = sigma[1]^2, (1, 2) = sigma[1]*rho[1, 2]*sigma[2], (1, 3) = sigma[1]*rho[1, 3]*sigma[3], (2, 1) = sigma[1]*rho[1, 2]*sigma[2], (2, 2) = sigma[2]^2, (2, 3) = sigma[2]*rho[2, 3]*sigma[3], (3, 1) = sigma[1]*rho[1, 3]*sigma[3], (3, 2) = sigma[2]*rho[2, 3]*sigma[3], (3, 3) = sigma[3]^2}))$

(7)

>

Download 3_gaussian.mw

Why not...

Answered: mmcdara 7891

April 05 2023

0 0

restart:
r := 4*t1+4*sqrt(t1^2-4*t2):
cond := t1 < 0:
piecewise(cond, S(cond), S(eval(cond, t1=-t1)))

                   /                   /      1   2\ \
          piecewise|t1 < 0, {t2 < 0}, { t2 <= - t1  }|
                   \                   \      4    / /

Technical problems and how to avoid them...

Answered: mmcdara 7891

April 03 2023

2 5

What technical problems?

Given the shape of the surface expr=f(k, r), the sulution found by NonlinearFiit is extremely doubtful.

You fill find a quite detailed analysis of your problem in the attached worksheet:

a recall of @tomleslie's solution,
a rapid evaluation of the stability of the corresponding best fit: "What does happen if the initial point is changed?",
an alternative approach based on Optimzation:-Minimize which,
1. when applied to the same expr, leads to a better agreement between the data and the firring curve,
2. when applied to the mean values (see details in the attached file), also leads to a close result which, at least theoritically is less sensitive to roundoff errors,
3. after having observed that the surface expr=f(k, r) contains 3 discontinuity locci, enables locating the valey where the minimum is to be seeked,
4. when applied over this valley and after proper re-parameterizations, gives an even deeper minimum.

fit_mmcdara.mw

And how to avoid them?

All the problems related to the discontinuities of the objective function can be removed with an ad hoc re-parameterization of the problem:

Let t the name of the regressor and y the name of the dependent variable.
Define tau as log(t)
Define z as 1/y

The only singularity is then located at k=0

With this double transformation, coupled (not necessary) with the "mean trick" one get a unique solution without any problem:

>

restart;

NATIVE DATA & MODEL

>

with(plots):
alias(ST = Statistics):

dat := Matrix([
               [1.0, 60], [1.0, 70], [1.0, 80], [1.0, 90],
               [1.0, 100], [1.0, 110], [1.0, 120], [1.0, 150],
               [2.0, 70], [2.0, 80], [2.0, 90], [2.0, 100], [2.0, 120],
               [2.0, 130], [2.0, 140], [2.0, 150], [2.0, 160], [2.0, 170],
               [3.0, 100], [3.0, 110], [3.0, 120], [3.0, 130], [3.0, 140], [3.0, 150],
               [3.0, 160], [3.0, 170], [3.0, 180], [3.0, 200], [3.0, 210], [3.0, 220]
             ]):

expr := k/(1 + (k/80 - 1)*2.73^(r*t)):

RE-PARAMETERIZATION OF THE ORIGINAL PROBLEM

>

# In order to reduce roundoff errors, it's better to replace all the couples by
# just 3 couples [1.0, m1], [2.0, m2], [3.0, m[3]] where each m[1], m[2] and m[3]
# represent the mean value of the data at locations 1.0, 2.0 &nd 3.0 .
#
# Doing this we lose some information and we have to define a new objective function
# J_mean which weights the 3 residuals of squares by the number of times each "t site"
# occurs.
#
# Theory says that the minima of J and J_mean, as their minimizers too, should be the same.
# But using J_mean instead of J will reduce roundoff errors:

L := convert(dat, listlist):
L_mean := [ seq([x, ST:-Mean(map2(op, 2, (select((u -> op(1, u)=x), L))))], x in {map2(op, 1, L)[]}) ];
W_mean := rhs~(ST:-Tally(map2(op, 1, L)));

(1)

>	y = k/(1 + (k/80 - 1)2.73^(rt)); eval(%, t=log(tau)); simplify(%) assuming r > 0; map((x -> 1/x), %);

y = k/(1+((1/80)*k-1)*2.73^(r*t))

y = k/(1+((1/80)*k-1)*2.73^(r*ln(tau)))

y = k/(1.+0.125e-1*tau^(1.004301609*r)*k-tau^(1.004301609*r))

1/y = (1.+0.125e-1*tau^(1.004301609*r)*k-tau^(1.004301609*r))/k

(2)

>	# Re-parameterization # tau = log(t) where represents the first column of dat # z = 1/x where represents the second column of dat L_mean_Re := map(x -> [exp(x[1]), 1/x[2]], L_mean)

(3)

>	# Define the model MODEL : tau ---> z model := 1/expr: model := simplify(eval(model, t=log(tau))): MODEL := unapply(model, tau): MODEL(tau);

(1+((1/80)*k-1)*tau^(1.004301609*r))/k

(4)

>

# The objective function J to minimize is the weighted residual sum of squares.

J_mean := ( add(W_mean *~ (MODEL~(map2(op, 1, L_mean_Re)) - map2(op, 2, L_mean_Re))^~2) );
J_mean := simplify(J_mean) assuming r > 0;

plot3d(J_mean, k=0..10, r=-1..1, style=surface, contours=50, axis[3]=[mode=log]);

opt_mean := Optimization:-Minimize(J_mean, { k >= 0, r >= 0}, initialpoint={k=1, r=0})

8*(-HFloat(0.010256410256410256)+(1+((1/80)*k-1)*2.718281828^(1.004301609*r))/k)^2+12*(-HFloat(0.008264462809917356)+(1+((1/80)*k-1)*7.389056099^(1.004301609*r))/k)^2+10*(-HFloat(0.006349206349206349)+(1+((1/80)*k-1)*20.08553692^(1.004301609*r))/k)^2

(5)

>	display( [ plot(eval(eval(expr, tau=exp(t)), opt_mean[2]), t=(min..max)(dat[..,1]), color=red), plot(L, style=point, symbol=solidcircle, symbolsize=16, color=blue) ] );

>	FIT := unapply( eval(eval(expr, tau=exp(t)), opt_mean[2]), t); ResidualSumOfSquares := add( (FIT~(dat[.., 1]) - dat[.., 2])^~2 );

(6)

>

Download fit_1_mmcdara.mw

MaxFlow doesn't seem useful...

Answered: mmcdara 7891

March 31 2023

1 5

You can get the answer using ShortestPath instead:

restart
with(GraphTheory):
with(SpecialGraphs):
g:=IcosahedronGraph():
HighlightVertex(g, {1,7}, green):
G := CopyGraph(g):
d := NULL:

while IsConnected(G) do
  s := ShortestPath(G, 1, 7);
  c := s[2..-2]:
       HighlightEdges(G, {seq({s[i], s[i+1]}, i=1..nops(s)-1)}, red);
  d := d, DrawGraph(G);
  G := DeleteVertex(G, c);
  print(s);
end do:
d := d, DrawGraph(G):

plots:-display(`<|>`(d))
                           [1, 2, 7]
                           [1, 3, 7]
                          [1, 5, 6, 7]
                          [1, 9, 8, 7]
                       [1, 4, 10, 12, 7]

ThisWay.mw

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

7891 Reputation

22 Badges

MaplePrimes Activity

These are answers submitted by mmcdara

You probably need to add assumptions...

Don't know if it may help ...

I don't have time right now to answe...

Rescaling is not the very problem...

This does not add much to what has been ...

A way...

For fun...maybe?...

plots:-intersectplot...

restart:with(Statistics):X := Vector([1,...

Why do you use a deprecated function?Not...

A lot of errors:...

How to correlate N gaussian random varia...

Why not...

Technical problems and how to avoid them...

MaxFlow doesn't seem useful...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again