mmcdara

3810 Reputation

17 Badges

6 years, 115 days

MaplePrimes Activity


These are replies submitted by mmcdara

@tomleslie 

Since Maple 18???
So I was as ignorant of this as Kitonum.

By the way, while reading again more carefully the post about the Maple's 2021 release I just saw that the new feature I talked about was titled  "Numerically Solve Vectorised Differential Equations".
My apologies

@vv 

There are even simpler examples where NextZero gives weird results:
 

# A correct result

f := x -> mul((x-1/k), k=1..4):
f(x);
s := RootFinding:-NextZero(f, 0);
for j from 2 to 4 do 
  s := RootFinding:-NextZero(f, s);
end do;
                          0.2499999999
                          0.3333333333
                          0.4999999999
                          0.9999999990

# And a stupid one

f := x -> mul((x-1/k), k=1..10):
f(x);
s := RootFinding:-NextZero(f, 0);
for j from 2 to 10 do 
  s := RootFinding:-NextZero(f, s);
end do;
                         0.09999999990
                         0.09999999990
                         0.09999999990
                         0.09999999990
                         0.09999999990
                         0.09999999990
                         0.09999999990
                         0.09999999990
                         0.09999999990
                         0.09999999990

Note that 

fsolve(f(x), x=0);

returns the 10 zeros.

And last but not least:

f := x -> (x-1/5)*(x-1/10):
s := RootFinding:-NextZero(f, 0);
s := RootFinding:-NextZero(f, s);
                         0.09999999990
                         0.09999999990

Probably something that can be fixed by an ad hoc tunning of NextZero options (as you said) but which seems very strange.

@Kitonum 

I guess @tomleslie is using Maple 2021 (see here 214351-Unveiling-Maple-2021)

@Carl Love 

Thanks, I take it as a compliment.

Do you remember that about a year ago I asked a question about "typing" random variables?

In one of your answers you detailed how to define an "RV" type (short for Random Variable) and the first idea I had preparing myself to answer @Mikhail Drugov was use the construction of this type so that a linear combination of RV quantities is an RV quantity too.
I hoped (naively) that Y :=1/2*X1 + 1/2*X2 being of type RV would then imply that Y inherits the properties (attibutes?) of X1 and X2.

Intuitively I imagined that a type was a kind of class that gathered objects (Maple's objects ?) sharing some common methods (attributes ?). In short something close to what exists in Python or Java for instance.
But this does not seem to be the way it is and I do not understand the connection between types and objects (if any).

Your paragraph concerning modules and objects talks about subtleties I don't completely understand.
All the more so I had always limited myself to this sentence from the help page (In Maple, objects are a special type of module) to avoid digging the topic. Are Maple's objects objects close to Python's or Java's?

@Mikhail Drugov 

Here are a few complements about HazardRate
HazardRate_Complement.mw

@Mikhail Drugov 

As far as I can remember I've pointed out problems in

  • sampling a Normal RV with the default method (that is without explicitely defining it) !!!
    Quantiles far from the mean by 4 standard deviations are over sampled. This is due to a poor choice of the default sampling method (Ziggurat method) instead of the classical Box-Muller one, for instance.
    The cure is to use method=envelope.
    The full thread is here A serious problem with Statistics:-Sample()  (sand15 is my professional login)
  • Related : 230618-Problem-In-Sampling-Gaussian-Random-Variables  (unanswered)
     
  • Cauchy's distribution sampling 226461-A-Remark-About-StatisticsSampling-Methodenvelope
     
  • Chi square Goodness of Fit test 227530-Some-Details-About-ChiSquareSuitableModelTest- unanswered
    The main problem is that the number of class is not defined according to the standard methods and that this can give poor results.
    I remember reporting a bug concerning this test: the number of degrees of freedom is not correctly calculated when parameters are estimated from the data.
     
  • I also reported several issues concerning discrete random variables (Binomial, NegativeBinomial, ...) but I shoould do a deeper research on this site to dig them out.
     
  • If you use the GammaDistribution be careful: shape and rate parameters are interverted with respect to the standard definition (shape first and rate next, the opposite in Maple)


I've probably discovered several other "little issues" (not reported), obtained disappointing performances (compared to R), blah-blah-blah.

But Statistics is also sometimes a rather stunning package as soon as you want to do symbolic computations.

So don't think that I don't recommend its use and see only flaws in it.
I use it practically every day and what bothers me the most is that it is too limited to do REAL statistics (here again compared to R). So I spend a lot of time implementing specific methods that I can find in R, and when I get tired of it, I send the data to R, process it, and return to Maple.

 

@Carl Love 

I didn't ask if if a linear combination where a and b are rationals would "work", just that it was a question we should reasonnably ask ourselves.

So: 

  • one may a priori think that if one uses rationals "Maple will correctly compute PDF(a*X1 + (1-a)*X2, t)" , as you say...
  • but it is worth checking

The verification is immediate, but the result is surprising :

restart:
with(plots):
with(Statistics):
X1 := RandomVariable(BetaDistribution(4, 4)):
X2 := RandomVariable(BetaDistribution(4, 4)):
Z := convert(1/91*X1 + 90/91*X2, rational):
plot(PDF(Z, t), t=0..1, view=[89/91..1, default]);

evalf(Variance(X1));
evalf(Variance(Z)); 


 

@acer 

Have you saw that writting W := value(PDF(0.75*X1 + 0.25*X2, t, inert=true)): gives a wrong result whatever the way you expand it, 
while W := value(PDF(3/4*X1 + 1/4*X2, t, inert=true)): without taking special precautions.
This is what I pointed out in my reply to Carl.
PDF.mw

IMO I would say there is some weakness in the way PDFs are computed for compound RVs.
I guess Maple uses a general strategy to dothat, but this trategy could be enhanced in the case a RV is a linear combination if two others (basically a convolution product of scaled pdfs).

The file below contains the numerical computation of the pdf of Z=a1*X1+a2*X2 in the case X1 and X2 are RVs whose pdfs are countinuous and have bounded supports.
There remain some integration problems when a1 (or a2) goes to 0 because the support of a1*X1 (a2*X2) goes to +oo, which implies a special care to the integration method.
Nevertheless it could be a starting point for a more reliable method than the one Maple uses, to compute the pdf of Z=a1*X1+a2*X2 ...

PDF_Linear_comb_rvs.mw

@Mikhail Drugov 

You can look my reply to @Carl Love

@Carl Love 

My first idea was to send a reply by saying that X1/2+ X2/2 worked well and that there was no need to increase Digits, just as you said.
But the things are not that simple...
As changing 0.5 to 1/2 "works", one can ask ourselves  if, for instance, a*X1 + (1-a)*X2 still "works" when a is a rational between 0 and 1 (infering this doesn't work if a is a floating point number between 0 and 1)?
Another question is: as the pdf of the sum of two independent RVs is equal to the convolution product of their pdfs, does replacing PDF (a*X1 + (1-a)*X2) by a convolution integral "works better"?

If you're interested to dig these questions deeper you will find some elements in the attached file
PDF.mw

@tomleslie 

Hi, 
This is indeed another way to solve the problem.
I'm going to look at it right now.
Thanks

SOLUTION FOUND!!! 



# bounded_sol denote the piecewise solution in my original post
#
# before using the code below you must unassigned s


sys := { 
  diff(v(t), t)=(s(t)-x(t)) * a(t), 
  diff(x(t), t)=v(t) * a(t), 
  s(t)=10*cos(t),
  beta_0(t)=piecewise(s(t)   >= 0, 1, 0),
  beta_2(t)=piecewise(s(t)-2 <= 0, 1, 0)-1,
  v(0)=0, 
  x(0)=0,
  a(0)=1
}:

evs  := [   
          [[x(t)-0, v(t) < 0], [a(t)=0, v(t)=0]],  
          [[x(t)-2, v(t) > 0], [a(t)=0, v(t)=0]],  
          [side(beta_0(t), pre), a(t)=1],
          [side(beta_2(t), pre), a(t)=1]
        ]:

esol := dsolve(sys, numeric, events=evs, discrete_variables=[a(t)], range=0..50):

display(
  odeplot(esol, [t, x(t)], 0..50, color=black, refine=1),
  bounded_sol, 
  gridlines=true, 
  size=[1200, 400]
)



The issue no longer exists.

@tomleslie 

No, it's not the same case.
In what you did, as soon as the ball touches the wall its velocity is reversed and the ball almost instantly moves backward.
It's easy a thing to model and simulates
What I want to simulate is completely different:

You are driving a Bugatti at low speed and suddenly you step on the gas.
You find yourself stuck against the seat and you stay there as long as your acceleration has not decreased enough: you don't bounce towards the windshield like the ball!
Or you push the brake pedal down and crash into the windshield (forget the seatbelt) and stay here until the deceleration drops or you hit the gas (miracles do happen sometimes).

 

More seriously, you have a closed cylindrical box with a ball inside, held in the center of the box by two springs of the same properties.
You want to trigger a signal as soon as the box is accelerated axially beyond a given value (whatever the direction) and maintain this signal as long as the acceleration has not fallen below this value.
A kind of inertial accelerometer.

@dharr 

Looks like a modern variation of Georgia O'Keefe's Flowers serie :-)

@Kitonum 

In the solution above, nowhere are the forces 70 lb and 60 lb taken into account....
... this wasn't asked, I just answered the question about angles.


Concerning your solution: I do not completely agree.

It seems that AC and AB are threads (with no flexion rigidity) and thus are pulls.
The support AO seems to be a rigid rod and without doing high mathematics it's clear this rod is a push.
The orientation of forces AC and AB means that yangtheary considers the material point A "isolated from the outside" and in equilibrium under the effect of 4 forces: thus the compression force in rod AO must be oriented from O to A.
For a simple planar problem the "triangle of forces" looks like this one

A solution where forces are only defined by their magnitudes and their lines of action, (AC, AB and OA lines) is not a complete solution.
For this to be the case, it is also necessary to provide the orientation of the force (if I had rotated the green arrow above by Pi the point A would no longer have been in equilibrium).
This is why I found that your solution is incomplete (using VectorCalculus would have done something better).

When you write "The sum of all forces applied to point  A  should equal 0." this is not "completely" true because this sentence has a meaning in terms of vectors only :"The sum of all forces applied to point  A  should be the null vector".

Last point: nothing is said concerning the type of joint between the rod AO and the wall at point O.
If it is a ball joint without any friction you can solve the problem and compute the forces.
If it is an embedding you have an hyperstatic problem that cannot be solved that easily and you must account for the flexion rigidity of the rod AO.




 

First 19 20 21 22 23 24 25 Last Page 21 of 87