sand15

1449 Reputation

15 Badges

11 years, 80 days

MaplePrimes Activity


These are answers submitted by sand15

UPDATED


The bug you observed is more subtle than what you think and does not reduce to a regression from one version to another.

Indeed, even if some Maple versions prior to Maple 2024 may return a no-error result (see @acer for instance), this does not guarantee this result is correct.
For instance Maple 2015 generates no error in calculating the pdf of sin(theta) but the resul it provides is wrong.
 

This latter claim is easy to diagnose:

  • If you know a few of Statistics it is obvious that both cos(theta) and  sin(theta) must have the same PDF... and clearly Maple does not provide identical results.
    Note that Maple 2015 computes correctly the pdf of cos(theta), but this is only by chance and Maple 2015 is not that lucky in computing the one of sin(theta)).
     
  • If you know a little bit more in Statistics, sin(theta) has, by definition, an arsine distribution , and if you know what thecorresponding pdf looks like, you can easily verify that Maple 2015's result is wrong.

The fact Maple 2015 (and likely all other Maple versions) returns (may return?) a wrong result is the consequence of a flawed mechanism for calculating PDF(𝜑(theta)) when 𝜑 is not a one-to-one map over the whole support of theta (see my comment ahead).
 

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

with(Statistics)

theta := RandomVariable(Uniform(0, 2*Pi))

_R

(2)

phi := arccos(-1+2*RandomVariable(Uniform(0, 1))); phi := arccos(RandomVariable(Uniform(-1, 1)))

arccos(-1+2*_R0)

 

arccos(_R1)

(3)

r := RandomVariable(Uniform(0, 1))^(1/3)

_R2^(1/3)

(4)

SinPhi := sin(phi)

(-_R1^2+1)^(1/2)

(5)

f__sp := unapply(PDF(SinPhi, t), t); f__sp(t)

piecewise(t <= 0, 0, 0 < t, 2*t*piecewise(t^2 < 0, 0, t^2 <= 1, piecewise(And(0 <= t^2, -t^4 < 0), piecewise(t^2 <= 0, 0, t^2 < 1, (1/4)/(-t^2+1)^(1/2), t^2 = 1, (1/4)/(t^2-1)^(1/2), 1 < t^2, 0), 0)+piecewise(And(0 <= t^2, -t^4 < 0), piecewise(t^2 < 0, 0, t^2 < 1, (1/4)/(-t^2+1)^(1/2), t^2 = 1, (1/4)/(t^2-1)^(1/2), 1 < t^2, 0), 0), 0))

(6)

SinTheta := sin(theta)

sin(_R)

(7)

f__st := unapply(PDF(SinTheta, t), t); f__st(t)

piecewise(t < 0, 0, t = 0, (3/2)/Pi, t < 1, 1/(Pi*(-t^2+1)^(1/2)), 1 <= t, 0)

(8)

CosPhi := cos(phi)

_R1

(9)

f__cp := unapply(PDF(CosPhi, t), t); f__cp(t)

piecewise(t < -1, 0, t < 1, 1/2, 0)

(10)

CosTheta := cos(theta)

cos(_R)

(11)

f__ct := unapply(PDF(CosTheta, t), t); f__ct(t)

piecewise(t <= -1, 0, t < 1, 1/(Pi*(-t^2+1)^(1/2)), 1 <= t, 0)

(12)

 

Supp := `~`[Support]([SinPhi, CosPhi, SinTheta, CosTheta], output = 'range')

[0 .. 1, -1 .. 1, 0 .. 1, -1 .. 1]

(13)

plots:-display(plot(f__sp(t), t = Supp[1], color = red, legend = typeset(sin('Phi'))), plot(f__cp(t), t = Supp[2], color = red, linestyle = 3, legend = typeset(cos('Phi'))), plot(f__st(t), t = Supp[3], color = blue, legend = typeset(sin('Theta'))), plot(f__ct(t), t = Supp[4], color = blue, linestyle = 3, legend = typeset(cos('Theta'))))

 

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
)

proc (t) options operator, arrow; piecewise(t <= -1, 0, t < 0, 2/(Pi*(-t^2+1)^(1/2)), 0 <= t, 0) end proc

 

proc (t) options operator, arrow; piecewise(t < 0, 0, t < 1, 2/(Pi*(-t^2+1)^(1/2)), 1 <= t, 0) end proc

 

 

# 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
)

fpos := proc (t) options operator, arrow; piecewise(t < 0, 0, t < 1, 2/(Pi*sqrt(-t^2+1)), 1 <= t, 0) end proc

 

(1/2)*piecewise(-t < 0, 0, -t < 1, 2/(Pi*(-t^2+1)^(1/2)), 1 <= -t, 0)+(1/2)*piecewise(t < 0, 0, t < 1, 2/(Pi*(-t^2+1)^(1/2)), 1 <= t, 0)

 

 


Why doesn't Maple compute correctly the PDF of  SinTheta?

d

# Let X some random variable and f : X --> Y a C1 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 fX(x)d(x), where fX denotes the PDF of X.
# This can be rewritten fX(x)d(x) = fX(x)/|f'(x)|d(y).
# So the PDF of Y is simply fY(y) = fX(x)/|f'(x)|,
#
# Introducing the reciprocal function y = f(-1) (which exists given the assumptions
# above) one has x=y(y) and thus
#         fY(y) = fX(y(y))|y'(y)|
#               = ( fX(y(y)) )'
#               = ( FX(y(y)) )'   (FX denotes the CDF of X)
#               = ( FY(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 WX of X into sub-intervals wn such that the restriction
#        of f has a unique inverse over each wn,
#    (2) compute the contributions Cn = (FX+y)'(y) over each wn 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]
)

[.4667653390, .5823642379]

 

[2.559228416, 2.674827315]

 

 

# 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);
 

piecewise(t <= -1, 0, t < 1, 1/(Pi*sqrt(-t^2+1)), 1 <= t, 0)

 

 

# 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);

piecewise(t <= 0, 0, t < 1, 2/(Pi*sqrt(-t^2+1)), 1 <= t, 0)

 

 

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

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

piecewise(t <= -1, 0, t < 1, 1/(Pi*sqrt(-t^2+1)), 1 <= t, 0)

 

 

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

 

 

Download Basics_Maple2015.mw

f(x) is singular at x=0 meaning its graph C is not continuous.
It's clear hat max(f(x) = +oo and min(f(x)) = -oo
So there is obviously no circle which answers your problem.

So i suppose that you saying "where the curve attains its maximum and minimum" has to be interpreted as "where diff(f(x), x) = 0 and diff(f(x), x, x) <> 0", am I right?
If it is so, because f(x) is antisymmetric if the point (a, f(a)), a being assumed > 0, is the point where f'(a)=0, then we have f'(-a)=0 too.
If there were a circle tangent to C at points (a, f(a)) and (-a, -f(a)), then, by symmetry condition, its center would br the point (0, 0).
So this circle could be tangenf to the curve C at points (a, f(a)) and (-a, -f(a)) if and only if a=0... which is impossible because f(0) is undefined.
So there is no circle solution of your problem.

@acer has provided the correct answer  but I am not pleased with it because it is a workaround to a Maple's weakness (unfortunately a workaround you must use if x and y have less "sympathetic" distrbutions).

As a rule I personnaly prefer to try and apply firstly a log transformation (logz = log(x) + log(y)) and compute next the pdf of z=exp(logz).
This trick does not work for any couple of distributions, but when it does, and it happens to be case here, it avoids using artificial (because a pdf is defined over the whole real line and there is no real justification to limit it to some interval) assumptions.

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)

piecewise(Zeta <= ln(2), 0, Zeta <= ln(3), -exp(Zeta)*(-Zeta+ln(2)), Zeta <= 2*ln(2), exp(Zeta)*(ln(3)-ln(2)), Zeta <= ln(2)+ln(3), exp(Zeta)*(-Zeta+ln(2)+ln(3)), ln(2)+ln(3) < Zeta, 0)

(1)

Z := exp(LZ):

PDF(Z, z)

piecewise(z <= 2, 0, z <= 3, ln(z)-ln(2), z <= 4, ln(3)-ln(2), z <= 6, ln(2)+ln(3)-ln(z), 6 < z, 0)

(2)
 

 

Download MultiplyPDFFunctions_sand15.mw


You can adjust "comfort ranges" as you please:

All_plots_Question_sand15.mw




Or, with adapted vertical scales

 

 

proposition.mw




You may be wondering why I did not plot arrows as you requested?
Here is my reply arrows.mw

Unless you make your own arrows by hand (which can be a waste of time), it is likely that neither plots nor plottools will give you satisfactory results (at least they don't for me)

You ask for a "mixture" of gaussian random variables and compute their sum instead, which is COMPLETELY DIFFERENT.

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

You seemed to agree with  @Carl Love when replying him and said "Mathematica formula for the variance is Var(Z)=w*a^2+(1-w)*b^2": in fact this formula is correct iif the expectations of both random variables are equal. and fron if not.
So it would be a mistake to take ir at face value.

Finally , for some reason I did not investigate deeply enough, the first part of @janhardo worksheet provide a wrong result (curiously identical to the so-called MMA's) while the second part where it defines mixtureDist is therefore incorrect.

A simple reasoning proves the so-called MMA result is wrong in case the tandom variable have different expectations:

Consider 2 gaussian random variables with respective means -A and +A and common standard deviations s << A.
Then a 50% mixture of both of themleads to a meand equal to 0 and a variance of the order than A^2.
As the formula  only accounts on standard deviation and not of expactions it is obviously false

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.

See here the correct answer:

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

with(plots):
with(Statistics):


What is a mixture of Gaussian Random Variables (GRV)?

Example

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]);

.5077706252*2^(1/2)*exp(-(9/2)*(x+1)^2)+.3385137501*2^(1/2)*exp(-(9/2)*(x-1)^2)

 

 


A mixture of GRVs IS NOT a sum of GRVs

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])

1.659700209*exp(-.3461538458-8.653846156*x^2-3.461538460*x)

 

 


Expectations

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

SGRV__mean := Mean(SGRV)

-.2000000000

 

-.2

(2)

Why are the expectations the same?

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');

proc (x) options operator, arrow; a*G1__pdf(x)+b*G2__pdf(x) end proc

 

a*(Int(x*G1__pdf(x), x = -infinity .. infinity))+b*(Int(x*G2__pdf(x), x = -infinity .. infinity))

 

a*(Int(x*G1__pdf(x), x = -infinity .. infinity)) = a*Mean(G1)

 

b*(Int(x*G2__pdf(x), x = -infinity .. infinity)) = b*Mean(G2)

 

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

 

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

(3)


Variances

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

GRV__var := Variance(SGRV);

1.071111111

 

0.5777777778e-1

(4)


Is MMA formula right (assuming you did not misinterpreted MMA's result...) ?

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

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

.1111111111

(5)


Observe that the variances of MGRV ans SGRV are different.
This is the result of the non linearity of the variance operator/

So, and this was already clear looking to the pdfs, a mixture is not a sum.

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]);

_m4479461408

(6)

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);

(3/2)*2^(1/2)*exp(-(9/2)*(x+1)^2)*w/Pi^(1/2)+(3/2)*2^(1/2)*exp(-(9/2)*(x-1)^2)*(1-w)/Pi^(1/2)

 

(1/2+(1/2)*erf((3/2)*(x+1)*2^(1/2)))*w+(1/2+(1/2)*erf((3/2)*(x-1)*2^(1/2)))*(1-w)

 

-2*w+1

 

4*w-4*w^2+1/9

 

-216*w*(2*w^2-3*w+1)/(-36*w^2+36*w+1)^(3/2)

 

-3*(1296*w^4-2592*w^3+1800*w^2-504*w-1)/(36*w^2-36*w-1)^2

(7)


Sampling

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

_m4510812096

(8)

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);

HFloat(-2.1959715574515415)

 

HFloat(2.1602521346759844)

 

 

[HFloat(-0.1996052599456992), HFloat(1.068640318864312)]

(9)


Mode

The mode can be quite difficult to catch given the wide variety of pdf a mixture can have.
So this is a quite simple, likely not optimal, procedure to do this.

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])

-.9999999797

(10)
 

 

Download Mixture_versus_Sum.mw

The correct mathematical form of an IC at infinity is limit(y(x), x=+oo) = 0 and the fact you can't write the IC this way in Maple does not authorize you to write y(+oo) = 0 instead.

The correct way to proceed is:

restart

kernelopts(version);

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

ode := diff(y(x),x$2)-2*diff(y(x),x)+y(x)=4*exp(-x);
IC  := y(+infinity)=0;
sol := dsolve([ode,IC])

diff(diff(y(x), x), x)-2*(diff(y(x), x))+y(x) = 4*exp(-x)

 

y(infinity) = 0

 

y(x) = -exp(I*Im(x))*signum(_C1)*infinity

(2)

IC  := y(L)=0;
sol := simplify(dsolve([ode,IC]))

y(L) = 0

 

y(x) = -L*exp(x)*_C1+exp(x)*x*_C1-exp(x-2*L)+exp(-x)

(3)

# The result below corresponds to what Maple did to get result (2)

limit(sol, L=+infinity)

y(x) = -exp(I*Im(x))*signum(_C1)*infinity

(4)

# Now consider result (3) and compute the limit this way

limit(limit(sol, x=L, left), L=+infinity)

limit(y(L), L = infinity) = 0

(5)

eval(sol, x=L);

y(L) = 0

 

y(L) = 0

(6)

# And below this is why you get an undefined result in your most recent question

eval(sol, L=+infinity);
eval(%, x=+infinity);

y(x) = -infinity*exp(x)*_C1+exp(x)*x*_C1-exp(x-infinity)+exp(-x)

 

y(infinity) = undefined

(7)
 

 

Download Limits.mw

Your "However, how I can create a cumulative histogram with corresponding polygon employing this same information?" question is unclear.
What do you mean by cumulative histogram?

Do you mean "an histogram of the cumulated data"... that is a collection of rectangles which would approximate the empirical cumulative distribution function (note I do not use the word "histogram" for this collection is not an histogram)?

If I interpret correctly what you seem to have in mind, here it is.
If I'm wrong, please clarify your question.

Done with Maple 2015 (I verified that Maple 18 contains all the features for file 1 to work with this version).

restart

with(Statistics):

Dados:=[-9.43, -4.42, -4.04, -2.88, -1.90, -1.81, -1.20, -1.16, -1.03, -.14, 1.27, 1.72, 1.97, 1.98, 2.24, 3.24, 3.64, 3.8, 5.1, 5.52]

[-9.43, -4.42, -4.04, -2.88, -1.90, -1.81, -1.20, -1.16, -1.03, -.14, 1.27, 1.72, 1.97, 1.98, 2.24, 3.24, 3.64, 3.8, 5.1, 5.52]

(1)

bounds := [-10, -7, -4, -1, 2, 5, 8]

[-10, -7, -4, -1, 2, 5, 8]

(2)

domains := [seq(bounds[i]..bounds[i+1], i=1..numelems(bounds)-1)]

[-10 .. -7, -7 .. -4, -4 .. -1, -1 .. 2, 2 .. 5, 5 .. 8]

(3)

T := TallyInto(Dados, domains, bins=1);

[HFloat(-10.0) .. HFloat(-7.0) = 1, HFloat(-7.0) .. HFloat(-4.0) = 2, HFloat(-4.0) .. HFloat(-1.0) = 6, HFloat(-1.0) .. HFloat(2.0) = 5, HFloat(2.0) .. HFloat(5.0) = 4, HFloat(5.0) .. HFloat(8.0) = 2]

(4)

X := [seq((bounds[i]+bounds[i+1])/2, i=1..numelems(bounds)-1)]

[-17/2, -11/2, -5/2, 1/2, 7/2, 13/2]

(5)

# It's better to use ListTools:-PartialSums than Statistics:-CumulativeSum because
# the latter returns HFloats numbers (unless you use 'UseHardwareFloats := false')
# that you will have to convert into integer for future use

Y := ListTools:-PartialSums(rhs~(T))

[1, 3, 9, 14, 18, 20]

(6)

WeightedX := [seq(X[i]$Y[i], i=1..numelems(X))]

[-17/2, -11/2, -11/2, -11/2, -5/2, -5/2, -5/2, -5/2, -5/2, -5/2, -5/2, -5/2, -5/2, 1/2, 1/2, 1/2, 1/2, 1/2, 1/2, 1/2, 1/2, 1/2, 1/2, 1/2, 1/2, 1/2, 1/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 7/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2, 13/2]

(7)

plots:-display(
  plot(
    X, Y
    , color=blue
    , thickness=3
  )
  , Histogram(
      WeightedX
      , binbounds = bounds
      , frequencyscale = absolute
      , color=gray
      , transparency=0.5
    )
  , view=[(min..max)(bounds), 0..max(Y)]
)

 
 

 

Cumulative_Histogram_with_Polygon_Line__Way_1.mw

In this Cumulative_Histogram_with_Polygon_Line_Way_2.mw worksheet the histogram is buit bars-by-bars using plottools:-rectangle instead of Statistics:-Histogram.
This avoids constructing the WeightedX list, which might be useful when Dados is a large sample, and offers more flexibility to produce smarter histigrams.

Illustration (note that the tickmarks could hve been displayed this way using the previous wotksheet).

So make your choice.

NOTE: the deprecated stats package (Maple 2015 and beyond) contained a useful feature enabling to draw histograms of weighted: Cumulative_Histogram_with_Polygon_Line__Way_1_stats_pckg.mw

As you use these two terms in the same question I understand you mean "envelope in the signalprocessing sense".

There exists a simple analogic device named Envelope Detector which, I quote Wikipedia,  "takes a (relatively) high-frequency signal as input and outputs the envelope of the original signal".

It is quite easy to draw the differential equation which models the input-output relation for this device (in the attached file Rd represents the resistance named R in the Wiki article, Cd the capacitance named C in this article and rd is an auxiliary resistance left to the diode and which Wiki does not account for).

This auxiliary edo is added to your eqn1 and eqn2 differential equations.
To manage the two "charge" and discharge" stages of the envelope detector device I use a discrete variable GT(t) whose value +1 corresponds to the charge of the cacitance and value -1 to its discharge.
The values (-1 and +1) are controlled by events in the dsolve function.

Because V2(t) is rapidly oscillating it is necessary to refine the odeplot outputs, which imposes in return to use the dsolve option range

File demodulation_toy_problem.mw shows how the envelope detector works. 
The solution depends on the three parameters Rd , Cd and rd which have to be adjusted. Each choice gives an envelope but some seem more "reasonable" than others, or at least seem to fit more than others a mind image of what an envelope should be.

File demodulation.mw applies this envelope detector to your problem.

You will likely observe a quite large display time (which goes with a large amount of memory used), but this is the price to pay to a fully resolved solution.
Which explains why I did not go up to t=6*T on my memory limited machine.

__________________________________________________________________________________


In those files the 3-parameter model  is a physical one.
A simpler 2-parameter physical model results in setting rd to some arbitrary value, for instance rd = 1 (in true envelope detector one has Rd >> rd ).
One may also think to a 2-parameter model (no longer that physical) whose parameters would be the two time-constants (one for the forward diode and the other for the reverse diode).

With this latter model the time-constant of the forward diode must be small enough for the envelope to follow faithfully the rise of V2(t).
If, as it seems it is your case, the successive maxima of V2(t) have increasing values (for t/T >~  0.9), then the time-constant for the reverse diode should be chosen large enough for  the envelope does not present significant decreasing between two consecutive peaks.

__________________________________________________________________________________


Finally,  demodulation.mw solves a 3-ode system but an alternative which could be interesting is to keep solving your initial 2-ode system first, and next dsolve the Envelope(t) ode alone while telling dsolve that V2(t) is now a known function (see known option and how to use it in dsolve/numeric help page).
This could mimic the post-processing stuff you mentioned.

Nevertheless not that concise as MMA

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

integrand:=sqrt(1+sin(x)^2);

(1+sin(x)^2)^(1/2)

(2)

with(IntegrationTools):

J0 := Int(sqrt(1+sin(x)^2), x)

Int((1+sin(x)^2)^(1/2), x)

(3)

J1 := Change(J0, sin(x)=u, u)

Int(-(u^2+1)^(1/2)*(-u^2+1)^(1/2)/(u^2-1), u)

(4)

J2 := value(J1);

EllipticE(u, I)

(5)

J3 := eval(J2, u=sin(x))

EllipticE(sin(x), I)

(6)

j3 := diff(J3, x)

cos(x)*(1+sin(x)^2)^(1/2)/(-sin(x)^2+1)^(1/2)

(7)

j4 := eval(j3, cos(x) = sqrt(1-sin(x)^2))

(1+sin(x)^2)^(1/2)

(8)

simplify(j4-integrand)

0

(9)
 

 

Download Maple_vs_MMA.mw

Try for instance

limit__FG1 := limit(FG1, epsilon = 0);

    exp(eta[1]) * exp(eta[2]) + exp(eta[1]) + exp(eta[2]) + 1

and have a look to Limit ( help(Limit) )

@C_R 

I was writting a quite detailed comment to your question but an unexpected box failure made me lost all of it.
So for short, because I'm going to bed right now, I send you a code with several rationl fraction fits of your data.
Jut tell me if you are interested in it and I will develop all of this tomorrow.

RationalFractionFit.mw

... and the answer remains the same: How can expect that solve will find the the solution of a non linear system of 17 equations in 8 unknowns?
Unless the rank of this system is 8 (and there is no guarantee solve will even succeed) his is a dead end.
So you will find in the attached file an alternative method which enables finding a solution "in the least squares sense" (something I already told you about a couple of times several months ago).

test-F-p_sand15.mw

As several members already told you the result is unevaluated in some "ancient" Maple versions.

Nevertheless evaluation can be obtained if you observe that the exponent 

2*sqrt(tau^2-1)/(sqrt(tau-1)*sqrt(tau+1))

is equal to 2 if tau > -1 and -2 otherwise and if you use proper simplification assumptions. 
Evaluation.mw

PS: I suppose someone smarter than me could achieve the same result in a simpler way.

try

Sol := unapply(sol, x):
eval(convert(IC, diff), Sol(x0));

Or in a row

eval(convert(IC, diff), unapply(sol, x)(x0))

Download A_way.mw

By the way, why are you so prudish that you refuse to call  this "other software" MMA ?

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