mmcdara

6279 Reputation

17 Badges

7 years, 325 days

MaplePrimes Activity


These are answers submitted by mmcdara



 

NULL

NULL

Animation Sinus

 

NULL

NULL

restart;

with(plottools):

with(plots):

h := 2:  

col := t -> `if`(floor(t/2/Pi)::even, red, blue);

proc (t) options operator, arrow; `if`((floor((1/2)*t/Pi))::even, red, blue) end proc

(1)

bg := display(
        circle([-2,0], 1, color=gray),
        pointplot([[-2,-h], [-2,h]], connect, color=black),
        pointplot([-2,0], symbol=solidcircle, symbolsize=10, color=blue),
        seq(
         plot(sin(t), t=(n-1)*2*Pi..n*2*Pi, color=col((n-1)*2*Pi), discont, gridlines)
         , n=1..4
       )
):

frame := proc(t)
        local P := [-2+cos(t), sin(t)], Q := [-2,sin(t)], H := [t,sin(t)];
        display(
                pointplot([[-2,0], P, Q], connect, color=blue),
                pointplot([P,Q,H], symbol=solidcircle, symbolsize=10, color=blue),
                line(Q, H, color=gray),
                arc([-2,0], 1, 0..t, color=col(t), thickness=5),
                line([-2,0], Q, color=col(t), thickness=5)
        );
end proc:

animate(frame, [t], t = 0 .. 8*Pi, background = bg, scaling = constrained, view = -h .. h, size = [1000, 700], frames = 120);

 

 

 

 


 

Download AnimSinusCouleurs_mmcdara.mw
 


For this kind of problems (or even more complex ones)  I like using the events option.
So I propose you two solutions:

  1. the first one uses only the events option to manage the changing of eq1's rhs,
  2. the second one uses options events and discrete_variables.
     

The first one is sufficient for your problem, but using also discrete_variables enables handling much more complex situations:

A_solution.mw

Of course this would have worked the same:

sol4 := dsolve(
  {
    diff(T(x), x) = piecewise(T(x) < 4, rhs1a, rhs1b),
    # or 
    # diff(T(x), x) = (rhs1b-rhs1a)*Heaviside(T(x)-4) + rhs1a,
    eqn2,
    T(0) = 0, 
    q(0) = 0
  },
  numeric
);

display(
  odeplot(sol , [x, T(x)], x = 1 .. 3.2, color=blue, legend="rhs1a"),
  odeplot(sol4, [x, T(x)], x = 1 .. 3.2, color=red , legend="rhs1a & rhs1b"),
  gridlines=true
);

but I prefer using events and discrete_variables which both offer more possibilities (which is of course a personal opinion).

Note that doing this would not give you the value of x where the "bifurcation" occurs (this is why the plotting range is set arbitrarily to 1..3.2). If this value of x maters to you, the simplest thing to do is to use events.

I'm not sure I inderstood correctly your problem, if not, just forget this.

restart

with(combinat):

# Holder

alph := parse~([$"A".."P"]);

# Number of players

NC   := 4:
game := alph[1..NC]:

# Random order of the players

parts := randperm(game):

# Build all series of NC/4 games each

game[1] := parts:
for n from 2 to NC-1 do
  game[n] := [game[n-1][1], ListTools:-Rotate(game[n-1][2..-1], 1)[]];
end do:

# Smart display

for __n from 1 to NC-1 do
  Games[__n] := op~([ seq([{game[__n][k..k+1][]}, vs,  {game[__n][k+2..k+3][]}], k=1..NC/4) ]);
end do;

[A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P]

 

[{A, B}, vs, {C, D}]

 

[{A, C}, vs, {B, D}]

 

[{A, D}, vs, {B, C}]

(1)

# Number of players

NC   := 16:
game := alph[1..NC]:

# Random order of the players

parts := randperm(game):

# Build all series of NC/4 games each

game[1] := parts:
for n from 2 to NC-1 do
  game[n] := [game[n-1][1], ListTools:-Rotate(game[n-1][2..-1], 1)[]];
end do:

# Smart display

for __n from 1 to NC-1 do
  Games[__n] := op~([ seq([{game[__n][k..k+1][]}, vs,  {game[__n][k+2..k+3][]}], k=1..NC/4) ]);
end do;

[{C, M}, vs, {F, L}, {C, L}, vs, {F, K}, {F, L}, vs, {J, K}, {F, K}, vs, {B, J}]

 

[{L, M}, vs, {F, K}, {F, L}, vs, {J, K}, {F, K}, vs, {B, J}, {J, K}, vs, {A, B}]

 

[{F, M}, vs, {J, K}, {F, K}, vs, {B, J}, {J, K}, vs, {A, B}, {B, J}, vs, {A, H}]

 

[{K, M}, vs, {B, J}, {J, K}, vs, {A, B}, {B, J}, vs, {A, H}, {A, B}, vs, {D, H}]

 

[{J, M}, vs, {A, B}, {B, J}, vs, {A, H}, {A, B}, vs, {D, H}, {A, H}, vs, {D, P}]

 

[{B, M}, vs, {A, H}, {A, B}, vs, {D, H}, {A, H}, vs, {D, P}, {D, H}, vs, {G, P}]

 

[{A, M}, vs, {D, H}, {A, H}, vs, {D, P}, {D, H}, vs, {G, P}, {D, P}, vs, {E, G}]

 

[{H, M}, vs, {D, P}, {D, H}, vs, {G, P}, {D, P}, vs, {E, G}, {G, P}, vs, {E, O}]

 

[{D, M}, vs, {G, P}, {D, P}, vs, {E, G}, {G, P}, vs, {E, O}, {E, G}, vs, {I, O}]

 

[{M, P}, vs, {E, G}, {G, P}, vs, {E, O}, {E, G}, vs, {I, O}, {E, O}, vs, {I, N}]

 

[{G, M}, vs, {E, O}, {E, G}, vs, {I, O}, {E, O}, vs, {I, N}, {I, O}, vs, {C, N}]

 

[{E, M}, vs, {I, O}, {E, O}, vs, {I, N}, {I, O}, vs, {C, N}, {I, N}, vs, {C, L}]

 

[{M, O}, vs, {I, N}, {I, O}, vs, {C, N}, {I, N}, vs, {C, L}, {C, N}, vs, {F, L}]

 

[{I, M}, vs, {C, N}, {I, N}, vs, {C, L}, {C, N}, vs, {F, L}, {C, L}, vs, {F, K}]

 

[{M, N}, vs, {C, L}, {C, N}, vs, {F, L}, {C, L}, vs, {F, K}, {F, L}, vs, {J, K}]

(2)
 

 

Download Games.mw

Are these lists of Games  ehaustive (because there are more with 16 players)?
 



EDITED 2023/03/19

Behind this intriguing title lies the fact that there exist several formulas to assess the weighted variance.


POINT 1
The most common formulas (biased and unbiased estimators) are given in the attached file (which is the one I already attached to my original answer)
Ecart_type.mw



POINT 2
This second file presents different formulas to assess the empirical weighted variance.
You will see that Maple's Variance(A, weights=B) corresponds to the formula one can find in two different external references.
Ecart_type.mw

I advice you to read carefully the two first pages of Berkeley where a discussion about what formula we should use is developped, in a nustshell this all depenf on the nature of the weights.
Given your weghts are integers, I would lean more towards saying that they represents the number of occurrences of each of your four outcomes?
Am I right?
If it so the good formula to use is the unbiased one in the first attached file here or the first one in the second attached file here. 
Whatever it would be great that Statistics:-Variance proposed different options to assess the variance of a weighted sample (just as it dioes to assess Quantiles)





POINT 3
In the previous files I didn't address the point concerning the discrepancy between Maple's Variance(A, weights=B) result and the result (variance=4.473676) your TI gives.

As you probably know these two points of view lead to different formulas for raw and centered moments of order larger than 1 (in particular the variance).
In this third file two points of view are discussed 

  1. In the first point of view the couple (A, B) is a sample of a random variable.
     
  2. In the second one (red text in the fle) the couple (A, B) represents the whole population.

At the end of it you will see that your TI adopts the second point of view.
Ecart_type.mw

 

See the edited answer for more details and explanations

 

 

restart:

with(Statistics):

A := [8, 9, 16, 18]:
B := [78, 31, 59, 81]:

# Direct formula fior the mean

Moyenne := add(A*~B) / add(B);
evalf(%);

3305/249

 

13.27309237

(1)

# Direct formula for the standard deviation
# (1): biased estimator

ET :=  sqrt(add((A -~ Moyenne)^~2*~B) / add(B));
evalf(%);

(1/249)*1240874^(1/2)

 

4.473675667

(2)

# Direct formula for the standard deviation
# (2): unbiased estimator

ET :=  sqrt(add((A -~ Moyenne)^~2 *~ B) / (add(B)-1));
evalf(%);

(1/15438)*4789153203^(1/2)

 

4.482686100

(3)

# Alternative way to compute the mean and the standard deviation (unbiased estimator)

C := [seq(A[i]$B[i], i=1..numelems(A))]:
[Mean, StandardDeviation](C);

[HFloat(13.273092369477911), HFloat(4.482686100195402)]

(4)
 

 

Download Ecart_type.mw
 

 

@Mariusz Iwaniuk @oggsait

... whose steps ar guided by the result the OP wants to get.
So it is not a solution that Maple produces by itself but a solution we force Maple to proivide...

In particular, step (5), I don't know how to force Maple to convert cosh(X)^2-1 into sinh(X)^2 without using eval (as I did) or simplify with the rule i use  in eval.
I also used several steps to transform sqrt(2)*sqrt(p[2])/sqrt(p[1]) into sqrt(2*p[2]/p[1]) because I wasn't able to do this more simply using ad hoc assumptions.

To @oggsait : I don't know if you will find some interest in that?
 

restart:

with(IntegrationTools):

f := 1/(u*sqrt(1+p*u^2/q))

1/(u*(1+p*u^2/q)^(1/2))

(1)

F := Int(f, u);

Int(1/(u*(1+p*u^2/q)^(1/2)), u)

(2)

`(1)`  = value(F);
`(2)`  = isolate(rhs(%)=X, u);
`(3)`  = allvalues(rhs(%));
`(4)`  = simplify([rhs(%)]) assuming p > 0, q < 0, x > x[0];
`(5)`  = eval(rhs(%), cosh = (z -> sqrt(1+sinh(z)^2)));
`(6)`  = eval(rhs(%), sinh(X)=S);
`(7)`  = simplify(rhs(%)) assuming S > 0;
`(8)`  = eval(rhs(%), S=sinh(X));
`(9)`  = eval(rhs(%), [p=P^2, q=Q^2]);
`(10)` = simplify(rhs(%)) assuming P > 0, Q < 0;
`(11)` = eval(rhs(%), Q=P*sqrt(q/p));
`(12)` = eval(rhs(%), sinh(X) = 1/sech(X));
`(13)` = eval(rhs(%), [X=x-x[0], q=``(-2*p[2]), p=p[1]]);  # assuming p[2] < 0

 

`(1)` = -arctanh(1/(1+p*u^2/q)^(1/2))

 

`(2)` = (u = RootOf(_Z^2*tanh(X)^2*p+tanh(X)^2*q-q))

 

`(3)` = (u = (-(tanh(X)^2*q-q)/(tanh(X)^2*p))^(1/2), u = -(-(tanh(X)^2*q-q)/(tanh(X)^2*p))^(1/2))

 

`(4)` = [u = (q/sinh(X)^2)^(1/2)/p^(1/2), u = -(q/sinh(X)^2)^(1/2)/p^(1/2)]

 

`(5)` = [u = (q/sinh(X)^2)^(1/2)/p^(1/2), u = -(q/sinh(X)^2)^(1/2)/p^(1/2)]

 

`(6)` = [u = (q/S^2)^(1/2)/p^(1/2), u = -(q/S^2)^(1/2)/p^(1/2)]

 

`(7)` = [u = q^(1/2)/(p^(1/2)*S), u = -q^(1/2)/(p^(1/2)*S)]

 

`(8)` = [u = q^(1/2)/(p^(1/2)*sinh(X)), u = -q^(1/2)/(p^(1/2)*sinh(X))]

 

`(9)` = [u = (Q^2)^(1/2)/((P^2)^(1/2)*sinh(X)), u = -(Q^2)^(1/2)/((P^2)^(1/2)*sinh(X))]

 

`(10)` = [u = -Q/(P*sinh(X)), u = Q/(P*sinh(X))]

 

`(11)` = [u = -(q/p)^(1/2)/sinh(X), u = (q/p)^(1/2)/sinh(X)]

 

`(12)` = [u = -(q/p)^(1/2)*sech(X), u = (q/p)^(1/2)*sech(X)]

 

`(13)` = [u = -(``(-2*p[2])/p[1])^(1/2)*sech(x-x[0]), u = (``(-2*p[2])/p[1])^(1/2)*sech(x-x[0])]

(3)

 


 

Download A_thirteen_steps_solution.mw


Even if it is possible to usre Runge-Kutta schems to solbe boundary value problems, this is not allowable in Maple.
So either you write your own RK algorithm or you trust Maple to use the scheme it founds most suited to.

solution.mw

The "ciulprit" is the term 

sqrt(gamma^2*sigma__d^2*sigma__e^4 + 4*sigma__e^2 + 4*sigma__v^2)

in the numerator of A.
The term under the radical is positive and can be noted T^2.
Now what is sqrt(T^2)?

b := sqrt(a^2):
simplify(b) assuming a > 0
                               a
simplify(b) assuming a < 0
                               -a

This is why you need assumptions.

Depending on the sign of T the numerator of A is equivalent to gamma^3 or gamma^5 when
gamma --> +oo (while the denominator of A is equivalent to gamma^4).

A detailed answer is developed in the attached file​​​​​​ ​limits_signum_explained.mw

To get less confusing notations a good practice is to assign capital letters to tandom variables.
So change delta by Delta, nu by Nu and so on. The relations will be clearer.

There is no need no write nu__1-nu__01  where nu_1 is a random variable and nu__01 its mean. The only advantage in doing this is to add an unnecesary  complexity!

Lastly, why not just ask Maple what those covariances are?

 

restart

with(Statistics):

RVS := [
         Nu__1    = RandomVariable(Normal(0, sigma__nu)),
         Nu__2    = RandomVariable(Normal(0, sigma__nu)),
         Delta__1 = RandomVariable(Normal(0, sigma__d)),
         Delta__2 = RandomVariable(Normal(0, sigma__d)),
         Delta__3 = RandomVariable(Normal(0, sigma__d3))
       ]

[Nu__1 = _R, Nu__2 = _R0, Delta__1 = _R1, Delta__2 = _R2, Delta__3 = _R3]

(1)

X__1 := beta__1*(Nu__1+Nu__2)+alpha__1*Delta__1+alpha__2s*Delta__2;
X__2 := beta__2*(Nu__1+Nu__2)+alpha__2*Delta__2+alpha__1s*Delta__1;
X__3 := beta__3*(Nu__1+Nu__2)+alpha__3*Delta__3;

beta__1*(Nu__1+Nu__2)+alpha__1*Delta__1+alpha__2s*Delta__2

 

beta__2*(Nu__1+Nu__2)+alpha__2*Delta__2+alpha__1s*Delta__1

 

beta__3*(Nu__1+Nu__2)+alpha__3*Delta__3

(2)

A := X__1*(-lambda__1*X__1-lambda__1*Delta__1+Nu__1);
B := X__2*(-lambda__2*X__2-lambda__2*Delta__2+Nu__2);
C := X__3*(-lambda__3*X__2-lambda__3*Delta__3+Nu__1+Nu__2);

(beta__1*(Nu__1+Nu__2)+alpha__1*Delta__1+alpha__2s*Delta__2)*(-(beta__1*(Nu__1+Nu__2)+alpha__1*Delta__1+alpha__2s*Delta__2)*lambda__1-lambda__1*Delta__1+Nu__1)

 

(beta__2*(Nu__1+Nu__2)+alpha__2*Delta__2+alpha__1s*Delta__1)*(-(beta__2*(Nu__1+Nu__2)+alpha__2*Delta__2+alpha__1s*Delta__1)*lambda__2-lambda__2*Delta__2+Nu__2)

 

(beta__3*(Nu__1+Nu__2)+alpha__3*Delta__3)*(-(beta__2*(Nu__1+Nu__2)+alpha__2*Delta__2+alpha__1s*Delta__1)*lambda__3-lambda__3*Delta__3+Nu__1+Nu__2)

(3)

`Cov(A, B)` := Covariance(eval(A, RVS), eval(B, RVS)):

`Cov(A, B)` := simplify(`Cov(A, B)`);

8*beta__1^2*beta__2^2*lambda__1*lambda__2*sigma__nu^4+8*beta__1*beta__2*alpha__1*alpha__1s*lambda__1*lambda__2*sigma__d^2*sigma__nu^2+8*beta__1*beta__2*alpha__2*alpha__2s*lambda__1*lambda__2*sigma__d^2*sigma__nu^2+2*alpha__1^2*alpha__1s^2*lambda__1*lambda__2*sigma__d^4+4*alpha__1*alpha__1s*alpha__2*alpha__2s*lambda__1*lambda__2*sigma__d^4+2*alpha__2^2*alpha__2s^2*lambda__1*lambda__2*sigma__d^4+4*beta__1*beta__2*alpha__1s*lambda__1*lambda__2*sigma__d^2*sigma__nu^2+4*beta__1*beta__2*alpha__2s*lambda__1*lambda__2*sigma__d^2*sigma__nu^2+2*alpha__1*alpha__1s^2*lambda__1*lambda__2*sigma__d^4+2*alpha__1*alpha__1s*alpha__2s*lambda__1*lambda__2*sigma__d^4+2*alpha__1s*alpha__2*alpha__2s*lambda__1*lambda__2*sigma__d^4+2*alpha__2*alpha__2s^2*lambda__1*lambda__2*sigma__d^4-4*beta__1^2*beta__2*lambda__1*sigma__nu^4-4*beta__1*beta__2^2*lambda__2*sigma__nu^4-2*beta__1*alpha__1*alpha__1s*lambda__1*sigma__d^2*sigma__nu^2-2*beta__1*alpha__2*alpha__2s*lambda__1*sigma__d^2*sigma__nu^2-2*beta__2*alpha__1*alpha__1s*lambda__2*sigma__d^2*sigma__nu^2-2*beta__2*alpha__2*alpha__2s*lambda__2*sigma__d^2*sigma__nu^2+alpha__1s*alpha__2s*lambda__1*lambda__2*sigma__d^4-beta__1*alpha__1s*lambda__1*sigma__d^2*sigma__nu^2-beta__2*alpha__2s*lambda__2*sigma__d^2*sigma__nu^2+beta__1*beta__2*sigma__nu^4

(4)

# and so on for Cov(A, C) and Cov(B, C)

 

 

Download No_problem.mw



One problem comes from the use of gamma, which is a protected name (Euler constant, see help page).
So when you write gamma, Maple will not consider it as an unknown in the solve process (which is why your solutions never contain something of the form gamma < (or >) something.

Another problem comes from using the assumptions directly in solve.
If you do this (look to my file for the notations), you will get solutions which contain b < 0

solve(expr > 0) assuming hyp:

# one non licit solution is 

{E = E, G = 1/(3*b), a = a, s = s, b < 0, -(-s+3*b)*E < c}

So you have to "filter" the solutions solve returns.

(
NOTES:

  1. I used  SOL_gt := solve(expr > 0) assuming hyp  but  SOL_gt := solve(expr > 0) seems to give the same result.
  2.  SOL_gt := solve(expr > 0, useassumptions=true) assuming hyp returns a wrong result (Maple 2015).
  3.  SOL_gt := solve({expr > 0, hyp}) seems to run endlessly  (Maple 2015).

)


Here is the solution:
 

restart

# gamma is the Euler constant

gamma;
evalf(gamma);

gamma

 

.5772156649

(1)

expr := (-(4*((E*s-2*c)*b+E*s*b))*b*G^2+(8*b^2*E+4*E*s*b+(4*(E*s-3*a*(1/2)-(1/2)*c))*b)*G-b*E-E*s+2*a-E*(b+s))/(2*(2*b*G-1)^2);

(1/2)*(-4*((E*s-2*c)*b+E*s*b)*b*G^2+(8*b^2*E+4*E*s*b+4*(E*s-(3/2)*a-(1/2)*c)*b)*G-b*E-E*s+2*a-E*(b+s))/(2*G*b-1)^2

(2)

hyp := a > 0, b > 0, c > 0, a > 2*c, G > 0, 2*b*G > a/c, s > 0, E >= 0

0 < a, 0 < b, 0 < c, 2*c < a, 0 < G, a/c < 2*b*G, 0 < s, 0 <= E

(3)

SOL_gt := solve(expr > 0) assuming hyp:
print~([%]):
 

{E = E, c = c, s = s, G < (1/2)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), b < 0}

 

{E = E, c = c, s = s, G < (1/3)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), b < 0, (1/2)/b < G}

 

{E = E, G = (1/3)/b, a = a, s = s, b < 0, -(-s+3*b)*E < c}

 

{E = E, c = c, s = s, b < 0, (1/3)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}

 

{E = E, G = G, b = 0, c = c, s = s, E*s < a}

 

{E = E, c = c, s = s, 0 < b, G < (1/3)/b, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}

 

{E = E, G = (1/3)/b, a = a, s = s, 0 < b, -(-s+3*b)*E < c}

 

{E = E, c = c, s = s, 0 < b, G < (1/2)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), (1/3)/b < G}

 

{E = E, c = c, s = s, 0 < b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), (1/2)/b < G}

(4)

# Eliminate solutions that do not verify somple hypothesis of the form
# "0 < X" ot "0 <= X"

simplehyp := select(has, {hyp}, 0);

good_gt := NULL:

for z in {SOL_gt} do
  map(x -> (coulditbe(x) assuming simplehyp[]), [z[]]);
  if not(member(false, %)) then
    good_gt := good_gt, z
  end if:
end do:

good_gt := {good_gt}:

numelems({SOL_gt}), numelems(good_gt);

print(cat('_'$140)):

print~(good_gt):

{0 <= E, 0 < G, 0 < a, 0 < b, 0 < c, 0 < s}

 

9, 4

 

____________________________________________________________________________________________________________________________________________

 

{E = E, G = (1/3)/b, a = a, s = s, 0 < b, -(-s+3*b)*E < c}

 

{E = E, c = c, s = s, 0 < b, G < (1/3)/b, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}

 

{E = E, c = c, s = s, 0 < b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), (1/2)/b < G}

 

{E = E, c = c, s = s, 0 < b, G < (1/2)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), (1/3)/b < G}

(5)

# I couldn't find a way to check if these four solutions verify the initial hypothesis (hyp)

SOL_lt := solve(expr < 0) assuming hyp:
print~([%]):
 

{E = E, c = c, s = s, G < (1/2)/b, b < 0, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}

 

{E = E, c = c, s = s, G < (1/3)/b, b < 0, (1/2)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}

 

{E = E, G = (1/3)/b, a = a, s = s, b < 0, c < -(-s+3*b)*E}

 

{E = E, c = c, s = s, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), b < 0, (1/3)/b < G}

 

{E = E, G = G, b = 0, c = c, s = s, a < E*s}

 

{E = E, c = c, s = s, 0 < b, G < (1/3)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1)}

 

{E = E, G = (1/3)/b, a = a, s = s, 0 < b, c < -(-s+3*b)*E}

 

{E = E, c = c, s = s, 0 < b, G < (1/2)/b, (1/3)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}

 

{E = E, c = c, s = s, 0 < b, (1/2)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}

(6)

# Eliminate solutions that do not verify somple hypothesis of the form
# "0 < X" ot "0 <= X"

good_lt := NULL:

for z in {SOL_lt} do
  map(x -> (coulditbe(x) assuming simplehyp[]), [z[]]);
  if not(member(false, %)) then
    good_lt := good_lt, z
  end if:
end do:

good_lt := {good_lt}:


print(cat('_'$140)):

print~(good_lt):

____________________________________________________________________________________________________________________________________________

 

{E = E, G = (1/3)/b, a = a, s = s, 0 < b, c < -(-s+3*b)*E}

 

{E = E, c = c, s = s, 0 < b, G < (1/3)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1)}

 

{E = E, c = c, s = s, 0 < b, (1/2)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}

 

{E = E, c = c, s = s, 0 < b, G < (1/2)/b, (1/3)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}

(7)

# check that no goo_gt solution is a good_lt solution and reciprovcally

good_gt union good_lt:
numelems(%);

good_gt intersect good_lt:
numelems(%);

8

 

0

(8)

``


 

Download Inequalities_mmcdara_1.mw

 


Is it too much to ask that you respect those who gave you the answer?

restart

alias(Q = Student:-Precalculus:-CompleteSquare):

M, P  := solve(_Z^2 + ``(epsilon - 1)*_Z - epsilon + theta, _Z):

op(1, M) + Q(expand(op(2, M)), epsilon);
op(1, P) + Q(expand(op(2, P)), epsilon);

-(1/2)*``(epsilon-1)+(1/2)*((epsilon+1)^2-4*theta)^(1/2)

 

-(1/2)*``(epsilon-1)-(1/2)*((epsilon+1)^2-4*theta)^(1/2)

(1)

Download A_way.mw


The first point is : "How do you define the mean and variance of g?
Look here reply_0.mw and read the last comment at the end of my answer.


Assuming that you want expressions for Ey and Vy which depend on the thetas, there is no need to define the latters as random variables before computing expressions Ey and Vy.
Finally your formula for Vy does not represents the variance of g but its 2nd order raw moment:
 

restart:

with(Statistics):


If you expect expressions of Ey and Vy to depend upon the thetas you must not define them before computing these expressions.
Instead of what you will get two numbers.

x[1] := RandomVariable(Normal(theta[1], theta[3]));

x[2] := RandomVariable(Normal(theta[2], theta[4]));

_R

 

_R0

(1)

g := 1 - (x[1] - 1)^2/9 - (x[2] - 1)^3/16;

1-(1/9)*(_R-1)^2-(1/16)*(_R0-1)^3

(2)

MyEy := CodeTools:-Usage( expand(Mean(g)) );
MyVy := CodeTools:-Usage( expand(Variance(g)) );

memory used=6.08MiB, alloc change=32.00MiB, cpu time=126.00ms, real time=1.26s, gc time=0ns

 

137/144-(1/9)*theta[3]^2-(1/9)*theta[1]^2+(2/9)*theta[1]-(1/16)*theta[2]^3-(3/16)*theta[4]^2*theta[2]+(3/16)*theta[2]^2+(3/16)*theta[4]^2-(3/16)*theta[2]

 

memory used=6.97MiB, alloc change=0 bytes, cpu time=78.00ms, real time=204.00ms, gc time=0ns

 

(9/64)*theta[4]^4+(2/81)*theta[3]^4+(9/64)*theta[2]^2*theta[4]^4-(9/32)*theta[2]*theta[4]^4+(4/81)*theta[1]^2*theta[3]^2-(8/81)*theta[1]*theta[3]^2+(4/81)*theta[3]^2+(15/256)*theta[4]^6+(9/256)*theta[4]^2*theta[2]^4-(9/64)*theta[4]^2*theta[2]^3+(27/128)*theta[4]^2*theta[2]^2-(9/64)*theta[4]^2*theta[2]+(9/256)*theta[4]^2

(3)

# Compare to what you want

Ey := -(1/9)*theta[1]^2 + (2/9)*theta[1] + 137/144 - (1/16)*theta[2]^3 + (3/16)* theta[2]^2 - (3/16)*theta[2] - (3/16*(-1+theta[2]))*theta[4]^2 - (1/9)*theta[3 ]^2:

Vy := (1/20736)*(9*theta[2]^3 + 16*theta[1]^2 - 27*theta[2]^2 - 32*theta[1] + 27*theta[2] - 137)^2 + (15/256)*theta[4]^6 + (45/256*(theta[2]^2 - 2*theta[2] + 1))*theta[4]^4 + ( 1/768*(45*theta[2]^4 + 32*theta[1]^2*theta[2] - 180*theta[2]^3 - 32*theta[1]^2 - 64*theta[ 1]*theta[2] + 270*theta[2]^2 + 64*theta[1] - 436*theta[2] + 301))*theta[4]^2 + (1/27)*theta[3]^4 + (1/216)*theta[3]^2*(3*theta[2]^3 + 16*theta[1]^2 - 9*theta[2]^2 - 32*theta[ 1] + 9*theta[2] - 35) + (1/24)*theta[3]^2*(-1+theta[2])*theta[4]^2;
 

(1/20736)*(9*theta[2]^3+16*theta[1]^2-27*theta[2]^2-32*theta[1]+27*theta[2]-137)^2+(15/256)*theta[4]^6+(45/256)*(theta[2]^2-2*theta[2]+1)*theta[4]^4+(1/768)*(45*theta[2]^4+32*theta[1]^2*theta[2]-180*theta[2]^3-32*theta[1]^2-64*theta[1]*theta[2]+270*theta[2]^2+64*theta[1]-436*theta[2]+301)*theta[4]^2+(1/27)*theta[3]^4+(1/216)*theta[3]^2*(3*theta[2]^3+16*theta[1]^2-9*theta[2]^2-32*theta[1]+9*theta[2]-35)+(1/24)*theta[3]^2*(-1+theta[2])*theta[4]^2

(4)

simplify(MyEy-Ey);
simplify(MyVy-Vy);

0

 

-(137/324)*theta[1]+(137/648)*theta[3]^2+(137/384)*theta[2]-(137/384)*theta[4]^2-18769/20736-(1/256)*theta[2]^6+(3/128)*theta[2]^5-(1/81)*theta[1]^4+(4/81)*theta[1]^3+(1/24)*theta[3]^2*theta[4]^2+(1/12)*theta[1]*theta[2]*theta[4]^2-(1/24)*theta[3]^2*theta[4]^2*theta[2]-(1/24)*theta[1]^2*theta[2]*theta[4]^2-(15/256)*theta[2]^4-(9/256)*theta[4]^4-(1/81)*theta[3]^4-(1/24)*theta[1]^2*theta[2]+(1/12)*theta[1]*theta[2]-(1/72)*theta[1]^2*theta[2]^3+(1/24)*theta[1]^2*theta[2]^2+(1/36)*theta[1]*theta[2]^3-(1/12)*theta[1]*theta[2]^2+(1/24)*theta[1]^2*theta[4]^2-(1/12)*theta[1]*theta[4]^2-(1/72)*theta[2]^3*theta[3]^2+(1/24)*theta[2]^2*theta[3]^2-(1/24)*theta[2]*theta[3]^2-(9/256)*theta[2]^2*theta[4]^4+(9/128)*theta[2]*theta[4]^4-(2/81)*theta[1]^2*theta[3]^2+(4/81)*theta[1]*theta[3]^2-(3/128)*theta[4]^2*theta[2]^4+(3/32)*theta[4]^2*theta[2]^3-(9/64)*theta[4]^2*theta[2]^2+(109/576)*theta[2]^3-(301/768)*theta[2]^2+(41/96)*theta[4]^2*theta[2]+(35/216)*theta[1]^2

(5)


You probably confounded the expression of the Variance with those of the Raw Moment of order 2:

M2 := Moment(g, 2);

(137/324)*theta[1]-(35/216)*theta[3]^2-(137/384)*theta[2]+(301/768)*theta[4]^2+18769/20736+(1/256)*theta[2]^6-(3/128)*theta[2]^5+(1/81)*theta[1]^4-(4/81)*theta[1]^3-(1/24)*theta[3]^2*theta[4]^2-(1/12)*theta[1]*theta[2]*theta[4]^2+(1/24)*theta[3]^2*theta[4]^2*theta[2]+(1/24)*theta[1]^2*theta[2]*theta[4]^2+(15/256)*theta[2]^4+(45/256)*theta[4]^4+(1/27)*theta[3]^4+(15/256)*theta[4]^6+(1/24)*theta[1]^2*theta[2]-(1/12)*theta[1]*theta[2]+(1/72)*theta[1]^2*theta[2]^3-(1/24)*theta[1]^2*theta[2]^2-(1/36)*theta[1]*theta[2]^3+(1/12)*theta[1]*theta[2]^2-(1/24)*theta[1]^2*theta[4]^2+(1/12)*theta[1]*theta[4]^2+(1/72)*theta[2]^3*theta[3]^2-(1/24)*theta[2]^2*theta[3]^2+(1/24)*theta[2]*theta[3]^2+(45/256)*theta[2]^2*theta[4]^4-(45/128)*theta[2]*theta[4]^4+(2/27)*theta[1]^2*theta[3]^2-(4/27)*theta[1]*theta[3]^2+(15/256)*theta[4]^2*theta[2]^4-(15/64)*theta[4]^2*theta[2]^3+(45/128)*theta[4]^2*theta[2]^2-(109/576)*theta[2]^3+(301/768)*theta[2]^2-(109/192)*theta[4]^2*theta[2]-(35/216)*theta[1]^2

(6)

simplify(M2-Vy);

0

(7)

VarianceFromM2 := M2 - MyEy^2:
simplify(VarianceFromM2 - MyVy)

0

(8)

 

 

Download reply_1.mw

 

 

Last comment: neither x[1] nor x[2] are random variables (in a mathematical sense) if the parameters of teir distributions are random variables.
But, for instance, x[1] conditionnally to  theta[1]=t[1] and theta[3]=t[3] is a gaussian random variables with parameters t[1] and t[3].
So do you want to compute conditional mean and variance or "whole" mean and variance?

A very likely non optimal code that could probably be enhanced.

restart

GeneratorsOf := proc(n)
  local N, F, G, s, t, r, d:

  uses combinat:

  N := 1:
  while fibonacci(N) < n do N := N+1 end do: N := max(6, N-1):

  F := {seq(fibonacci(i), i=2..N)};
  G := convert(convert~(powerset(F) minus {{}}, list), list):
  s := select((g -> add(g)=n), G);

  if numelems(s) = 1 then
    return s[1]
  else
    t, r := selectremove((x -> numelems(x) = 1), s);
    if t <> [] then
      return t[1]
    else
     d := map(e -> e[2..-1] -~ e[1..-2], r):
     zip((x, y) -> if not(has(y, 1)) then x end if, r, d);
     return select(has, %, {F[]})[1]
    end if:
  end if;

end proc:
  
  

GeneratorsOf(10);

[2, 8]

(1)

GeneratorsOf(100)

[3, 8, 89]

(2)

GeneratorsOf(142)

[1, 5, 13, 34, 89]

(3)

for n from 2 to 100 do
  printf("%4d  %a\n", n, GeneratorsOf(n))
end do;

   2  [2]
   3  [3]

   4  [1, 3]
   5  [5]
   6  [1, 5]
   7  [2, 5]
   8  [8]
   9  [1, 8]
  10  [2, 8]
  11  [3, 8]
  12  [1, 3, 8]
  13  [5, 8]
  14  [1, 13]
  15  [2, 13]
  16  [3, 13]
  17  [1, 3, 13]
  18  [5, 13]
  19  [1, 5, 13]
  20  [2, 5, 13]
  21  [8, 13]
  22  [1, 21]
  23  [2, 21]
  24  [3, 21]
  25  [1, 3, 21]
  26  [5, 21]
  27  [1, 5, 21]
  28  [2, 5, 21]
  29  [8, 21]
  30  [1, 8, 21]
  31  [2, 8, 21]
  32  [3, 8, 21]
  33  [1, 3, 8, 21]
  34  [13, 21]
  35  [1, 34]
  36  [2, 34]
  37  [3, 34]
  38  [1, 3, 34]
  39  [5, 34]
  40  [1, 5, 34]
  41  [2, 5, 34]
  42  [8, 34]
  43  [1, 8, 34]
  44  [2, 8, 34]
  45  [3, 8, 34]
  46  [1, 3, 8, 34]
  47  [13, 34]
  48  [1, 13, 34]
  49  [2, 13, 34]
  50  [3, 13, 34]
  51  [1, 3, 13, 34]
  52  [5, 13, 34]
  53  [1, 5, 13, 34]
  54  [2, 5, 13, 34]
  55  [21, 34]
  56  [1, 55]
  57  [2, 55]
  58  [3, 55]
  59  [1, 3, 55]
  60  [5, 55]
  61  [1, 5, 55]
  62  [2, 5, 55]
  63  [8, 55]
  64  [1, 8, 55]
  65  [2, 8, 55]
  66  [3, 8, 55]
  67  [1, 3, 8, 55]
  68  [13, 55]
  69  [1, 13, 55]
  70  [2, 13, 55]
  71  [3, 13, 55]
  72  [1, 3, 13, 55]
  73  [5, 13, 55]
  74  [1, 5, 13, 55]
  75  [2, 5, 13, 55]
  76  [21, 55]
  77  [1, 21, 55]
  78  [2, 21, 55]
  79  [3, 21, 55]
  80  [1, 3, 21, 55]
  81  [5, 21, 55]
  82  [1, 5, 21, 55]
  83  [2, 5, 21, 55]
  84  [8, 21, 55]
  85  [1, 8, 21, 55]
  86  [2, 8, 21, 55]

  87  [3, 8, 21, 55]
  88  [1, 3, 8, 21, 55]
  89  [34, 55]
  90  [1, 89]
  91  [2, 89]
  92  [3, 89]
  93  [1, 3, 89]
  94  [5, 89]
  95  [1, 5, 89]
  96  [2, 5, 89]
  97  [8, 89]
  98  [1, 8, 89]
  99  [2, 8, 89]
 100  [3, 8, 89]

 
 

 

Download fib.mw

In red in the attached file

TriagleSemblablesTest.mw
 

An1 := animate(display, [('f')(T2, phi)], phi = 0 .. arctan(4/3), background = T1);
...
An2 := animate(display, [('g')(T3, h)], h = 0 .. -15.8, background = T1);

 

1 2 3 4 5 6 7 Last Page 2 of 53