# Your problem will be simpler to solve if you consider Z0 instead of G0

# convert(Z01, exp) replaces hyperbolic trigs by exponentials

T := simplify(convert(Z01, exp));

# replace T = A^N/2 * B^N/2 by (A*B)^N/2
# and discard the exponent already in common with G0

T0 := op([1, 1], T)* op([2, 1], T)

# the question becomes: Is T0 equal to Z0?

simplify(T0 - Z0)

restart

kernelopts(version)

branches(arccot)

# The inverse function of 'cot' depends on the branch you consider

`#msup(mo("cot"),mo("-1"))` := solve(cot(y)=A, y, allsolutions);

about(_Z1);

Originally _Z1, renamed _Z1~:
  is assumed to be: integer

# Main branch

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=0)

# branch (-Pi, 0)

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=-1)

# branch (-2*Pi, -Pi)

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=-2)

# branch (Pi, 2*Pi)

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=1)

restart

# MathML hexa code of symbol x

`#mo("≃")`

# Use this symbol as a binary operator.. version 1

`≃`(Pi, 3.14);  # not very pretty

# Version 2;: define operator 'ae' (almost equal)
# see acer's comment
`print/ae` := proc(x,y)
   uses Typesetting;
   mrow(Typeset(x), mo("⁢"),
        mo("≃",fontweight="bold"), mo("⁢"),
        Typeset(y));
end proc:

ae(Pi, 3.14);  # prettier

`print/ra` := proc(x,y)
   uses Typesetting;
   mrow(Typeset(x), mo("⁢"),
        mo("→",fontweight="bold"), mo("⁢"),
        Typeset(y));
end proc:

plot(cos(x), x=0..2, title = typeset(ra(cos(x)=0, ae(x, 1.6))), titlefont=[Times, bold, 16])

ex:=(a+b*c)^2

shorten := table([_Inert_POWER="`^`", _Inert_SUM="`+`", _Inert_PROD="`*`", "("="[", ")"="]"]):

inertex := ToInert(ex):

shortstr := convert( eval(inertex, {_Inert_NAME = (u -> parse(u)), _Inert_INTPOS = (u -> u)}), string):

shortstr := copy(shortstr):
for i in [indices(shorten, nolist)] do
  shortstr := StringTools:-SubstituteAll(shortstr, i, short[i])
end do:

cat("[", shortstr, "]")

restart

kernelopts(version)

DocumentTools:-Tabulate(
  [
    plots:-display(
      Histogram(Sample(SinPhi  , 10^4), title=typeset(sin('Phi')), style='polygon', transparency=0.8),
      plot(f__sp(t), t = Supp[1], color = red)
    ),
    plots:-display(
      Histogram(Sample(CosPhi  , 10^4), title=typeset(cos('Phi')), style='polygon', transparency=0.8),    
      plot(`f__cp`(t), t=Supp[2], color=red)
    ),
    plots:-display(
      Histogram(Sample(SinTheta, 10^4), title=typeset(sin('Theta')), style='polygon', transparency=0.8),    
      plot(`f__st`(t), t=Supp[3], color=red),
      plots:-textplot([0, 10, "WRONG PDF"], font=[Tahoma, bold, 20], color=red)
    ),
    plots:-display(
      Histogram(Sample(CosTheta, 10^4), title=typeset(cos('Theta')), style='polygon', transparency=0.8),    
      plot(`f__ct`(t), t=Supp[4], color=red)
    )
  ]
)

# Correct computation of PDF(SinTheta)

Tneg := RandomVariable(Uniform(-Pi, 0)):
Tpos := RandomVariable(Uniform(0, Pi)):

fneg := unapply(PDF(sin(Tneg), t), t);
fpos := unapply(PDF(sin(Tpos), t), t);

plots:-display(
  plot((1/2)*~[fneg, fpos](t), t=-1..1, color=[red, blue]),
  Histogram(Sample(SinTheta, 10^5), title=typeset(sin('Phi')), style='polygon', transparency=0.8, minbins=100)
  , gridlines=true
)

# Or, directly from symmetry considerations

Tpos := RandomVariable(Uniform(0, Pi)):
fpos := unapply(PDF(sin(Tpos), t), t);

f__SinTheta := 1/2 * (fpos(-t) + fpos(t));

plots:-display(
  plot(f__SinTheta, t=-1..1, color=[red, blue]),
  Histogram(Sample(SinTheta, 10^5), title=typeset(sin('Phi')), style='polygon', transparency=0.8, minbins=100)
  , gridlines=true
)

d

# Let X some random variable and f : X --> Y a C¹ function of X.
# Let Supp(X) the support of X.
#
# In the simplest case where phi has a unique inverse for each x ∊ Supp(X), then
# an infinitesimal neihghborhood V(x) of diameter d(x) of x has the image
# d(y) = |f'(x)|d(x) (d(x) assumed > 0 and f'(x) denotes the derivative of f(x)).
#
# The probability that an outcome of X belongs to the interval [x-d(x), x+d(x)]
# has the value f_X(x)d(x), where f_X denotes the PDF of X.
# This can be rewritten f_X(x)d(x) = f_X(x)/|f'(x)|d(y).
# So the PDF of Y is simply f_Y(y) = f_X(x)/|f'(x)|,
#
# Introducing the reciprocal function y = f^(-1) (which exists given the assumptions
# above) one has x=y(y) and thus
#         f_Y(y) = f_X(y(y))|y'(y)|
#               = ( f_X(y(y)) )'
#               = ( F_X(y(y)) )'   (F_X denotes the CDF of X)
#               = ( F_Y(y) )'      (= derivative of the CDF of Y = f(x))
#
#
#
# For your problem it is worth noting that neither the sine nor cosine f functions have
# a unique inverse on the support of tandom variable theta.
#
# When f does not have a unique inverse one must proceed more carefully:
#    (1) split the whole support W_X of X into sub-intervals w_n such that the restriction
#        of f has a unique inverse over each w_n,
#    (2) compute the contributions C_n = (F_X+y)'(y) over each w_n and sum them to get the
#        final result.
#
# Here is a graphical illustration for the sine function.
# The probability that Y belongs to the interval [0.45, 0.55] is equal to the probability
# that X belongs to interval I1 OR that X belongs to interval I2.


I1 := sort([fsolve(sin(x)=0.45, x=0..Pi/2.), fsolve(sin(x)=0.55, x=0..Pi/2.)]);
I2 := sort([fsolve(sin(x)=0.45, x=Pi/2...Pi), fsolve(sin(x)=0.55, x=Pi/2...Pi)]);


plots:-display(
  plot(sin(x), x=0..2*Pi, color=blue, labels=["X", "Y"], labelfont=[Tahoma, bold, 12])
  , plot(0.5, x=0..2*Pi, color=black)
  , plot(0.55, x=0..2*Pi, color=black, linestyle=3)
  , plot(0.45, x=0..2*Pi, color=black, linestyle=3)
  , plottools:-rectangle([-0.05, 0.45], [0.05, 0.55], color=black)
  , plots:-textplot([2*Pi, 0.55, y+dy], align={left, above})
  , plots:-textplot([2*Pi, 0.45, y-dy], align={left, below})
  , plot([[I1[1], 0.45], [I1[1], 0]], color=black, linestyle=3)
  , plot([[I1[2], 0.55], [I1[2], 0]], color=black, linestyle=3)
  , plottools:-rectangle([I1[1], -0.05], [I1[2], 0.05], color=black)
  , plot([[I2[1], 0.45], [I2[1], 0]], color=black, linestyle=3)
  , plot([[I2[2], 0.55], [I2[2], 0]], color=black, linestyle=3)
  , plottools:-rectangle([I2[1], -0.05], [I2[2], 0.05], color=black)
  , plots:-textplot([I1[1], 0.1, x-dx], align={left})
  , plots:-textplot([I1[2], 0.1, x+dx], align={right})
  , plots:-textplot([I2[1], 0.1, x-dx], align={left})
  , plots:-textplot([I2[2], 0.1, x+dx], align={right})
  , size=[800, 400]
)

# A domain where phi (aka 'sin') has a unique inverse:

T := RandomVariable(Uniform(-Pi/2, Pi/2)):

PDF(sin(T), t);
plot(%, t=-1..1);
 

# But on the same domain phi (aka 'cos') does not have a unique inverse... so a wrong result

PDF(cos(T), t);
plot(%, t=-1..1);

# To get a correct result one must translate T in such a way that T_translated has a unique
# 'cos' inverse once the translation is done.
# Here, simply:

PDF(cos(T + Pi/2), t);
plot(%, t=-1..1);

# It should be clear to you that cos(theta) and sin(theta) (once computed correctly)
# must have the same distribution.

restart:

with(Statistics):

x := RandomVariable(Uniform(1, 2)):
y := RandomVariable(Uniform(2, 3)):

LX := log(x):
LY := log(y):

LZ := LX + LY:

PDF(LZ, zeta)

Z := exp(LZ):

PDF(Z, z)

restart

kernelopts(version)

with(plots):
with(Statistics):

G1 := RandomVariable(Normal(-1, 1/3)):
G2 := RandomVariable(Normal(+1, 1/3)):

p := 0.6:


prop := [p, 1-p]:

MGRV__pdf := add(PDF(G||i, x)*prop[i], i=1..2);

plot(MGRV__pdf, x=-3..3, title="Mixture of 2 GRVs\n", titlefont=[Verdana, bold, 12]);

SGRV := add(G||i*prop[i], i=1..2):

SGRV__pdf := PDF(SGRV, x);

plot(SGRV__pdf, x=-3..3, title="Sum of 2 GRVs\n", titlefont=[Verdana, bold, 12])

MGRV__mean := int(x*MGRV__pdf, x=-infinity..+infinity);

SGRV__mean := Mean(SGRV)

mgrf__pdf := x -> a * G1__pdf(x) + b * G2__pdf(x);

J0 := Int(x * mgrf__pdf(x), x=-infinity..+infinity):

with(IntegrationTools):

J1 := Expand(J0);


# obviously

op(1, J1) = a * :-Mean('G1');
op(2, J1) = b * :-Mean('G2');

# and thus

:-Mean(Mixture('G1', 'G2', 'p', 1-'p')) = 'p' * :-Mean('G1') + (1-'p') * :-Mean('G2');

# which, given the linearity of the :-Mean operator gives

:-Mean(Mixture('G1', 'G2', 'p', 1-'p')) = :-Mean('p' * 'G1' + (1-'p')*'G2');

MGRV__var := int(x^2*MGRV__pdf, x=-infinity..+infinity) - MGRV__mean^2;

GRV__var := Variance(SGRV);

MMA__var := p*Variance(G||1) + (1-p)*Variance(G||2);

# So the so-called MMA-formula you provide is incorrect.

GaussianMixture := proc(GRV::list(RandomVariable), W:=list(nonnegative))

  local N, V, pdf, cdf, m:

  N := numelems(GRV):

  if nops(W) <> N then
    error "GRV and P must have the same number of elements"
  end if:

  V := map(rv -> op(0, (attributes(rv)[3]):-ParentName()), GRV):
  if `not`(`and`(op(is~(V =~ NormalDistribution)))) then
    error "not all the GRV are Gaussian Random Variables"
  end if:
  
  pdf := unapply( add(PDF(GRV[i], x)*W[i], i=1..N), x);
  cdf := unapply( add(CDF(GRV[i], x)*W[i], i=1..N), x);
  m   := add(Mean(GRV[i])*W[i], i=1..N);


  Distribution(
    PDF = unapply(pdf(x), x)
    , CDF = unapply(cdf(x), x)

    # It is not mandatory to define these two latter quantities.
    # It remains up to you to do it or not.

    , Mean = m

    , Variance = int((x-m)^2 * pdf(x), x=-infinity..+infinity)
  )


end proc:
 

M := GaussianMixture([G1, G2], [w, 1-w]);

PDF(M, x);
CDF(M, x);
Mean(M);
Variance(M);

# Quantities not defined in the definition of a GaussianMixture random variable
# still remain computable.

Skewness(M);
Kurtosis(M);

Mn := GaussianMixture([G1, G2], [p, 1-p]);

QL := Quantile(Mn, 1e-4);
QU := Quantile(Mn, 1-1e-4);

sample := Sample(Mn, 10^4, method=[envelope, range=QL..QU]):

display(
  Histogram(sample, color="LightGray", style='polygon')
  , plot(PDF(Mn, x), x=(min..max)(sample), color=red, thickness=2)
);

[Mean, Variance](sample);

FindMode := proc(Mn)

  local QL, QU, f, q, h, sig, dec, dom, inc, maximizers, i, ix:

  QL := Quantile(Mn, 1e-4);
  QU := Quantile(Mn, 1-1e-4);

  f := unapply( diff(PDF(Mn, x), x), x):
  q := [seq(QL..QU, (QU-QL)/100)]:
  h := evalf(f~(q));

  sig := signum~(h):

  # Search for points corresponding to a change of sign in f(x)
  # For each sign changing at location 'loc', define the smallest neighbourhood
  # ('dom') within which f(x) gets both negative and positive values.
  dec := ListTools:-Search(-1, sig);
  dom := [dec-1, dec];

  inc := ListTools:-Search(+1, subsop(seq(i=0, i=1..dec-1), sig));
  while inc > 0 do
    dom := dom, [inc, inc+1];
    dec := ListTools:-Search(-1, subsop(seq(i=0, i=1..inc-1), sig));
    dom := dom, [dec-1, dec];
    inc := ListTools:-Search(+1, subsop(seq(i=0, i=1..dec-1), sig));
  end do:
  dom := [dom];

  # Solve f(x)=0 within each 'dom'ain of odd rank
  maximizers := NULL;
  for i from 1 to numelems(dom) by 2 do
    maximizers := maximizers, fsolve(f(x), x=q[dom[i][1]]..q[dom[i][2]])
  end do;
  maximizers := [maximizers]:
 
  # I assume below there is only one single mode (not to be confounded to maxima whose
  # number can be equal to the number of mixed random variables).
  # Accounting for the multiple mode situation only require to add a test in the structure
  # below.
  if numelems(dom) > 1 then
    ix := sort(map(i -> evalf(eval(PDF(Mn, x), x=i)), maximizers), output=permutation, `>`):
    return maximizers[ix[1]]
  else
    return maximizers[1]
  end if:

end proc:

FindMode(Mn, [G1, G2])