4019 Reputation

17 Badges

6 years, 177 days

MaplePrimes Activity

These are replies submitted by mmcdara

@Carl Love

I'm quite sure you smiled a lot when you got this plot !
As I'm sure you immediately saw where my mistake was and savored in advance the idea that I would be perplexed by this graph ?
After 20 seconds of total confusion, I realized that the plot range for p  was completely wrong.

Here is the corrected file.
Could you please run it with your recent Maple ?

PS: I do appreciate your joke

@Preben Alsholm 

Thank you Preben.
Glad to read I can obtain the closed form with Maple 2021.
I'll wait until next Monday when I'm in the office to use your tip (until then I'm at home with Maple 2015).

If it is not abusing your kindness, could you tell me if this works correctly with Maple 2021 or 2022?
(h and p should superimpose; p is the PDF of the product of V1 by V2)



U__1 := RandomVariable('Beta'(1, 1)):
U__2 := RandomVariable('Beta'(5, 5)):

V__1 := U__1 * 2 + 1:
V__2 := U__2 / StandardDeviation(U__2) / 3 + 3:

s__1 := Sample(V__1, 10^6):
s__2 := Sample(V__2, 10^6):
s    := s__1 *~ s__2:
h    := FrequencyPlot(s, minbins=1000, color=blue, transparency=0.4, title=typeset(pdf(''V__1''*''V__2''))):

support__1 := Support(V__1, 'output=range'):

f__1 := unapply(PDF(V__1, t), t):
f__2 := unapply(PDF(V__2, t), t):
f__W := z -> int(f__1(t)*f__2(z/t)*(1/abs(t)), t=support__1):


p := plot(simplify(f__W(z)), z=0..1, color=red):

display(p, h, size=[1000, 400])



Thanks in advance.
If you have no time to do that, just drop it. It will wait until Monday, the future of mankind is not at stake.


In case, like me, you are an airhead and forget to follow up on your questions, I left a comment on your last Post.
Sorry for having interfered here

In Maple 2015:

  • Select  View > Show/Hide contents
  • Switch off "Output"

Nothing alike in newer versions ?
Seems strange, look here

I agree that converting into Heaviside generates points at which the function is undefined, which in turns will appear as discontinuity points according to standard definition of continuity based on the use of limits (left_limit(f(x), x=c) = right_limit(f(x), x=c) = f(c)... which needs f(c) to be defined).

But your answer is not satisfying: If the conversion into Heaviside functions was the only reason which explains the occurrence of discontinuity points, then this would be the case whether the supports of the URVs are abstract or numeric. 
But as we can see at the bottom of this file:

  • the PDF of the sum of 2 abstract URV contains 3 discontinuity points,
  • but these same points are removed once the parameters of the URVs are set to numerical values.  

If converting into Heaviside functions before computing the convolution product xas the sole explanation of the apperance of discontinuity points (which I understand), why instanciating the parametres of the URVs to numerical values removes these discontinuity points ?



U1 := RandomVariable(Uniform(-1, 1)):
U2 := RandomVariable(Uniform(-1, 1)):

PDF(U1+U2, t);
discont(%, t)

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




# Using convolution
# (we get the same results as before)

f  := unapply(PDF(U1, t), t);
ff := int(f(z)*f(t-z), z=-infinity..+infinity);

discont(ff, z);

f := proc (t) options operator, arrow; piecewise(t < -1, 0, t < 1, 1/2, 0) end proc


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




# Now let's work with abstract uniform RV.
# Three discontinuity points appear

U1 := RandomVariable(Uniform(a, b)):
U2 := RandomVariable(Uniform(a, b)):

fa  := unapply(PDF(U1, t), t):
ffa := int(fa(z)*fa(t-z), z=-infinity..+infinity) assuming b > a:
ffa := simplify(ffa);

discont(ffa, t)

ffa := piecewise(a < b, piecewise(a < t-b, piecewise(b < t-b, 0, (2*b-t)/(-b+a)^2), piecewise(a < t-a, -(2*a-t)/(-b+a)^2, 0)), 0)


{2*a, 2*b, b+a}


# Once instanciated with a=-1 and b=+1 the discontinuity points are removed

g := simplify(eval(ffa, [a=-1, b=1]));
discont(g, t);

g := piecewise(t <= -2, 0, t <= 0, 1/2+(1/4)*t, t <= 2, 1/2-(1/4)*t, 2 < t, 0)




# What if PDF(U1, t) is converted into Heaviside functions?

U1 := RandomVariable(Uniform(a, b)):
U2 := RandomVariable(Uniform(a, b)):

fa  := unapply(convert(PDF(U1, t), Heaviside), t):
ffh := int(fa(z)*fa(t-z), z=-infinity..+infinity) assuming b > a:
ffh := simplify(ffh);

discont(ffh, t);   # 3 discontinuity points

# Once instanciated with a=-1 and b=+1 the discontinuity points are removed

h := simplify(eval(ffh, [a=-1, b=1]));
discont(h, t);   # discontinuity points removed

ffh := piecewise(b < a, 0, piecewise(a < t-b, piecewise(b < t-b, 0, (2*b-t)/(-b+a)^2), piecewise(a < t-a, -(2*a-t)/(-b+a)^2, 0)))


{2*a, 2*b, b+a}


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







@Joe Riel   I think you will better understand what I am doing if I attach an example.
Here is (not up to date) the worsheet dedicated to the building of a package which contains the construction of a few designs of experiments

I hope this will answer some of your questions.
I stay at your disposal for any further clarification.

@Joe Riel  @acer

you're partially right: I'm writing the code in a worksheet and used that method to break the code into chunks instead of one giant input region (structured programming). 
But I also write archives which gather a collection of procedures corresponding to a same domain (for instance a StatisticsTools package which contains procedures to build many design of experiments and others to realize specific data analysis).

To both of you, about the "style" I use to write modules,
When I started the project to write a huge Maple application to model, simulate and analyse the behaviour of complex mechanical systems, I realized that I had to do something else than just writing a  chunk of code in a worksheet.
Having learnt about modules, packages and Maple Librarie Archives, I posted a question on Mapleprimes to ask advices about the way to use them.
(I think this question, probably sent under my profesional login sand15  five ot six years ago, could be recover on this site).

By the time I used to do something inspired by what I found in the module help page:

z5 := module()
    export plus, times;
    plus := (a,b) -> a + b mod 5;
    times := (a,b) -> a * b mod 5;
end module;

If I'm not mistaken I've been advised to proceed differently:

  • either by doing what I'm doing now (procedures defined outside the module),
  • or by using $include (which looks like what we do in Fortran?)

I opted for the first way

@acer  @tomleslie

Thank you both for your answers.

To @acer:  what kind of troubles did I refer to?
The module M is built this way:

P1 := proc(a, b) a+b end proc:
P2 := proc(a, b) a*b + P1(a, b) end proc: 
P3 := proc(a, b) a*b + M:-P1(a, b) end proc:
M := module()
    option package:
    export P1 := :-P1,
           P2 := :-P2,
           P3 := :-P3;
end module:
M:-P2(3, 4), M:-P3(3, 4), 
                               19, 19

Here P2 and P3 give the same result.
But in some circumstances it was not the case and I needed to refer explicitly to the module name (that is write somethong like M:-Pn instead of Pn)
It's difficult for me to recover the exact situations when I needed to do that for this concerns the job I do at the office (which I cannot access from home). So I'll be browsing my archives next Monday to identify these situations and I  will keep you informed.

But I seem to remember that these situations appeared when a procedure was like this one

Pm := proc... end proc:

Pn := proc(...)
  local maplet:
  uses Maplets, Maplets:-Elements:
  maplet := Maplet(
  result := Display(maplet)
end proc

and that Ineeded to write 



But I don't want to waste any more of your time and I'll get back to you next Monday or Tuesday (if I find the damn problems again).


@Carl Love 

You're right Carl.
But I also use Maple 2021 at work and I'm aware of this build-in-interpreter. I even did little tests whit it.
OK, I can't say if this interpreter enables using specific Python packages (Maple 2015, sobs), but for what I remember I don't think so.

I'm aware that saying that it's mere hot air

@Carl Love 

You are right, my main concern was to explain @rockyicer why its code did not worked as expected.
Of course this can be donne in a much simpler way. I was about to publish something when I saw that you had taken the lead (no doubt in a more concise way).

I vote up

@Carl Love 

Great solution, I vote up.
I was completely unaware of the existence of this package. Nevertheless I'm not sure I quite understantthe meaning of the different points.


Thanks Christopher.

My first idea was indeed to use op.
I had a little bit more complex expression of `*` type, something like this

a := Term1*Term2*...*Unit(some_unit);
b := `*`(op(1..-2, a));  # expected to be the dimensionless version of a

Unfortunately I can't put the hand of the expression I was on, but I remember that Unit was not the rightmost operator and thus

`*`(op(1..-2, a))

didn't give me the dimensionless version of a.

It would be easier for anyone to help you if you provided your Maple code. 
Click on the big green color to load your code.

@Carl Love 

Thank you Carl for your involvement.


Thanks acer, this works perfectly well.
The case of b is indeed interesting, even if I never add meters and milimeters (but I guess some people could do so).

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