Advanced Search

mmcdara

7891 Reputation

22 Badges

9 years, 54 days
Contact mmcdara

MaplePrimes Activity


These are questions asked by mmcdara

Behaviour of print~...

mmcdara 7891 Maple 2015
print
February 26 2019
2 2

Hi, 

Could you explain me where the third output comes from?

print~([a, b]);
   a
   b
   [ ]

Thanks in advance

Inverse of the CDF of a Gamma RV...

mmcdara 7891 Maple
solve
February 20 2019
2 1

Hi, 

The CDF of a continuous random variable of support S is a bijective function f : S --> [0, 1].
So I expected that Maple would return only one solution when the command solve(f(t)=u, t) assuming u >=0, u <=1

When f(t) is the CDF of a Gamma random variable withparameters (2, 2), Maple returns two different solutions.

Could you explain me where the "spurious" solution (red curve) comes from?


 

> 

restart
> 

with(Statistics):
> 

C := RandomVariable(GammaDistribution(2, 2))

_R

 (1)
> 

f := unapply(CDF(C, t), t) assuming t > 0;

plot(f(t), t=0..10);

proc (t) options operator, arrow; 1-(1/2)*t*exp(-(1/2)*t)-exp(-(1/2)*t) end proc

  

 
> 

# f being a bijection from [0, +infinity) to [0, 1], its inverce does exist

s := solve({f(t)=u, t >=0}, t) assuming u >= 0, u <= 1;

{t = -2*LambertW((-1+u)*exp(-1))-2}, {t = -2*LambertW(-1, (-1+u)*exp(-1))-2}

 (2)
> 

pc := plot(CDF(C, t), t=0..10, color=blue):
T  := plottools:-transform((x, y) -> [y, x]):
plots:-display(T(pc), plot(rhs(op(s[1])), u=0..0.95, color=red), plot(rhs(op(s[2])), u=0..0.95, color=cyan, transparency=0.5, thickness=5));

 
> 

 


 

Download ICDF_of_GammDistribution.mw

How to display an Array filled with 0?...

mmcdara 7891 Maple
array lprint
February 16 2019
1 3

Hi, 

I'm stuck on this problem to which a careful reading of the help pages would probably give an answer:

Why is the return of Array(-1..1, [1$3]) of a different type than the one of Array(-1..1, [0$3]) ?

> 

restart:
> 

B := Array(-1..1,[1$3]);
lprint(B);

B := Array(-1..1, {(1) = 1})

  

Array(-1 .. 1, {-1 = 1, 0 = 1, 1 = 1})

  
> 

B := Array(-1..1,[0$3]);
lprint(B)

B := Array(-1..1, {(1) = 0})

  

Array(-1 .. 1, {})

  
> 

seq(B[n], n=-1..1)

0, 0, 0

 (1)
> 

 


 

Download Zero_Array.mw

A remark about Statistics:-Sampling(...,...

mmcdara 7891 Maple
sample statistics
February 15 2019
0 8

Hi, 

Recently a few questions concerning the sampling of the Cauchy distribution and the sampling of a truncated Normal distribution have been posted (mainly by  @jalale).
This post is concerned by the sampling of a truncated (standard) Cauchy distribution.

In a first part the efficiency fo two methods is adressed in the case of a non-truncated Cauchy distribution:

  • The "standard" Maple's command Statistics:-Sample(Cauchy(0, 1), N)
  • And a general method a priori very efficient if one knows the ICDF (Inverse Cumulative Function Distribution). It happens that this ICDF is just cot(U*Pi)  where U is a Uniform RV over [0, 1].
     

The second part adresses the sampling of a truncatedCauchy distribution with two methods:

  • The "standard" Maple's command Statistics:-Sample(Cauchy(0, 1), N, method=[envelope, range=...])
  • The method based on the use of the ICDF

 

Results:

Test1 (non-truncated Cauchy distribution) 

  • "Standard" Maples sampling outperforms the ICDF based method in terms of :
    • memory occupation: ICDF is twice more demanding
    • cpu time: ICDF is ten times slower
       

Test2 (truncated Cauchy distribution) 

  • ICDF based method i outperforms "Standard" Maples sampling oin terms of :
    • memory occupation: Maples "envelope sampling" method is twice more demanding
    • cpu time: Maples "envelope sampling" method is two times slower


But, beyond these simple observations, a disturbing problem is: the "envelope sampling" method seems to not return the correct distribution (at least when used this waymethod=[envelope, range=a..b]  with a < b)
This is confirmed by the two last plot where histogram and PDF are uperimposed.

Do you think this problem can be avoided by another parameterization of the "envelope sampling" method or that it reveals some underlying problem with it?

PS: I did not investigate further for other distributions .

 


 

> 

 

Sampling the Cauchy distribution

Maple's default sampling method outperformes the adhoc method

> 

restart
> 

with(Statistics):
> 

C := RandomVariable(Cauchy(0, 1))

_R

 (1)
> 

f := unapply(CDF(C, t), t);

proc (t) options operator, arrow; 1/2+arctan(t)/Pi end proc

 (2)
> 

finv := unapply(-solve(f(t)=u, t), u)

proc (u) options operator, arrow; cot(u*Pi) end proc

 (3)
> 

# "natural" way to proceed

N  := 10^6:

S1 := CodeTools:-Usage(Sample(C, N)):

memory used=7.71MiB, alloc change=39.63MiB, cpu time=69.00ms, real time=69.00ms, gc time=8.72ms

  
> 

# Let's try the sampling strategy based on the inverse of the CDF
# Usually it starts from sampling a Uniform RV on [0, 1] and
# next applies finv to the result.
#
# Smart but inefficient

U  := RandomVariable(Uniform(0., 1)):
S2 := CodeTools:-Usage(finv~(Sample(U, N))):

memory used=145.02MiB, alloc change=7.63MiB, cpu time=4.94s, real time=3.04s, gc time=2.61s

  
> 

# Much more efficient
#
# Given the special form of finv it's simpler to sample a Unirorm RV
# on [0, Pi] and apply "cot" to the result

pi := evalf(Pi):
U  := RandomVariable(Uniform(0., pi)):
S2 := CodeTools:-Usage(cot~(Sample(U, N))):

memory used=15.28MiB, alloc change=0 bytes, cpu time=652.00ms, real time=237.00ms, gc time=577.78ms

  

Sampling a truncated Cauchy distribution

Example 1:
with(Statistics) + method=[envelope, range=-10..10]

The adhoc method outperforms Maple's default sampling method

> 

S1 := CodeTools:-Usage(Sample(C, N, method=[envelope, range=-10..10])):

Histogram(S1);

memory used=8.88MiB, alloc change=-7.63MiB, cpu time=322.00ms, real time=260.00ms, gc time=93.77ms

  

 
> 

p  := Probability(C < -10, numeric);
q  := 1-Probability(C > +10, numeric);
U  := RandomVariable(Uniform(p*pi, q*pi)):
S2 := CodeTools:-Usage(cot~(Sample(U, N))):

Histogram(S2);

HFloat(0.03172551743055352)

  

HFloat(0.9682744825694465)

  

memory used=15.28MiB, alloc change=7.63MiB, cpu time=170.00ms, real time=113.00ms, gc time=92.11ms

  

 

Sampling a truncated Cauchy distribution

Example 2:
with(Statistics) + method=[envelope, range=-1..1]

The adhoc method outperformes Maple's default sampling method

> 

S1 := CodeTools:-Usage(Sample(C, N, method=[envelope, range=-1..1])):


scaling := Probability(C < +1, numeric) - Probability(C < -1, numeric);
plots:-display( Histogram(S1), plot(PDF(C, t)/scaling, t=-1..1, thickness=3, color=red) );

memory used=8.08MiB, alloc change=0 bytes, cpu time=246.00ms, real time=215.00ms, gc time=46.62ms

  

HFloat(0.5)

  

 
> 

p  := Probability(C < -1, numeric);
q  := 1-Probability(C > +1, numeric);
U  := RandomVariable(Uniform(p*pi, q*pi)):
S2 := CodeTools:-Usage(cot~(Sample(U, N))):

plots:-display( Histogram(S2), plot(PDF(C, t)/scaling, t=-1..1, thickness=3, color=red) );

HFloat(0.25)

  

HFloat(0.75)

  

memory used=15.28MiB, alloc change=7.63MiB, cpu time=163.00ms, real time=106.00ms, gc time=85.93ms

  

 
> 

 


 

Download CAUCHY_adhoc-sampling.mw

Maybe not the right place to tell about ...

mmcdara 7891 Maple
pde pdsolve differential-equations
February 02 2019
1 3

I'm presently interested in PDE and I have just discovered the impressive work Nasser Abbasi,has done, and keeps doing, concerning the solution of PDE benchmarks with Maple and Mathematica.
https://www.12000.org/my_notes/pde_in_CAS/pdse3.htm

It seems that the result Maple returns for the test case below is not the general solution.

  • 4.19 first order PDE of three unknowns
    problem number 19
    (from example 3.5.4, p 212, nonlinear ode’s by Lokenath Debnath, 3rd edition)


It is rather simple to see that any spherical function u(x, y, z) =f(x^2+y^2+z^2) is a solution of the PDE.
Then any function of the form u(x, y, z) = f(x^2+y^2+z^2) * exp(x+y+z) *C1 (C1 being any constant) is also a solution.
Maple returns only the solution u(x, y, z) = exp(C2*(x^2+y^2+z^2)) * exp(x+y+z) * C1


 

> 

restart:
> 

u := (x, y, z) -> f(x^2+y^2+z^2)

proc (x, y, z) options operator, arrow; f(x^2+y^2+z^2) end proc

 (1)
> 

expr := (y-z)*diff(u(x, y, z), x)+(z-x)*diff(u(x, y, z), y)+(x-y)*diff(u(x, y, z), z)

2*(y-z)*(D(f))(x^2+y^2+z^2)*x+2*(z-x)*(D(f))(x^2+y^2+z^2)*y+2*(x-y)*(D(f))(x^2+y^2+z^2)*z

 (2)
> 

simplify(%);

0

 (3)
> 

u := (x, y, z) -> C*f(x^2+y^2+z^2)*exp(x+y+z)

proc (x, y, z) options operator, arrow; C*f(x^2+y^2+z^2)*exp(x+y+z) end proc

 (4)
> 

expr := (y-z)*diff(u(x, y, z), x)+(z-x)*diff(u(x, y, z), y)+(x-y)*diff(u(x, y, z), z)

(y-z)*(2*C*(D(f))(x^2+y^2+z^2)*x*exp(x+y+z)+C*f(x^2+y^2+z^2)*exp(x+y+z))+(z-x)*(2*C*(D(f))(x^2+y^2+z^2)*y*exp(x+y+z)+C*f(x^2+y^2+z^2)*exp(x+y+z))+(x-y)*(2*C*(D(f))(x^2+y^2+z^2)*z*exp(x+y+z)+C*f(x^2+y^2+z^2)*exp(x+y+z))

 (5)
> 

simplify(%);

0

 (6)
> 

pdsolve((y-z)*diff(U(x, y, z), x)+(z-x)*diff(U(x, y, z), y)+(x-y)*diff(U(x, y, z), z)=0, U(x,y,z),'build')

U(x, y, z) = exp((1/2)*x^2*_C2)*exp(_C1*x)*exp((1/2)*y^2*_C2)*exp(_C1*y)*_C3*_C5*_C4*exp((1/2)*z^2*_C2)*exp(_C1*z)

 (7)
> 

combine(rhs(%), exp)

exp((1/2)*x^2*_C2+_C1*x+(1/2)*y^2*_C2+_C1*y+(1/2)*z^2*_C2+_C1*z)*_C3*_C5*_C4

 (8)
> 

 


 

Download Problem_19.mw

First 40 41 42 43 44 45 46 Page 42 of 48