mmcdara

772 Reputation

12 Badges

3 years, 106 days

MaplePrimes Activity


These are questions asked by mmcdara

Hi, 

Here's a problem that's been bothering me for a long time:

The Statistics[Distribution] help page describes the way to construct one's own Random Variable (RV).
All the examples construct a distribution foo := Distribution(...) next used to define the RV X by typing
X := RandomVariable(foo)

On the other way, constructing a RV Y from predefined Maple's distributions writes 
Y := RandomVariable(bar(sequence of parameters))

The "user defined RV feature" (MyRV := RandomVariable(foo)) seems to be restricted to distributions where all the parameters are numerically instanciated, contrary to "Maple's defined distributions" where these parameters can be kept formal).
So my question: Is it possible to define your own Distribution (module foo) in such a way that the definition of the RV X writes
X := RandomVariable(foo(sequence of parameters))

The attached file contains an example for the Levy distribution (the second part is just a workaround where a substitution is used to balance my inability to define a RV the way above)


 

restart:

with(Statistics):

A particular Levy distribution

Levy Distribution  Levy(m=0, s=1)
m : location parameter
s := dispersion parameter

m := 0:
s := 1:

pdf := unapply(piecewise(x<0, 0, sqrt(s/2/Pi)*exp(-s/2/(x-m))/(x-m)^(3/2)), x);
cdf := unapply(erfc(sqrt(s/2/(x-m))), x);


samplingMethod := proc(N)
                     m +~ s /~ ( Quantile~(Normal(0, 1), Sample(Uniform(1/2, 1), N)) )^~2
                  end proc;

proc (x) options operator, arrow; piecewise(x < 0, 0, (1/2)*2^(1/2)*exp(-(1/2)/x)/(Pi^(1/2)*x^(3/2))) end proc

 

proc (x) options operator, arrow; erfc((1/2)*2^(1/2)*(1/x)^(1/2)) end proc

 

proc (N) `~`[:-`+`](m, ` $`, `~`[:-`/`](s, ` $`, `~`[:-`^`](`~`[Statistics:-Quantile](Normal(0, 1), Statistics:-Sample(Uniform(1/2, 1), N)), ` $`, 2))) end proc

(1)

Levy := Distribution(
                      PDF=(x->pdf(x)),
                      CDF= (x->cdf(x)),
                      RandomSample = (N->samplingMethod(N))
                    );

_m4561537536

(2)

X := RandomVariable(Levy):
PDF(X, x);
Sample(X, 5);

piecewise(x < 0, 0, (1/2)*sqrt(2)*exp(-1/(2*x))/(sqrt(Pi)*x^(3/2)))

 

Vector[row]([.569870365734461, .357020836498979, 39.1445163385190, .340597592385078, 1.23209867072410])

(3)


How should I define a formal Levy distribution in such a way that, as for any predefined RV in Maple,
a particular instance is defined by  X := RandomVariable(Levy(M, S))  
where M and S are numerical values?

m := 'm':
s := 's':

pdf := unapply(piecewise(x<0, 0, sqrt(s/2/Pi)*exp(-s/2/(x-m))/(x-m)^(3/2)), x);
cdf := unapply(erfc(sqrt(s/2/(x-m))), x);
samplingMethod := proc(N)
                     m +~ s /~ ( Quantile~(Normal(0, 1), Sample(Uniform(1/2, 1), N)) )^~2
                  end proc;

proc (x) options operator, arrow; piecewise(x < 0, 0, (1/2)*2^(1/2)*(s/Pi)^(1/2)*exp(-(1/2)*s/(x-m))/(x-m)^(3/2)) end proc

 

proc (x) options operator, arrow; erfc((1/2)*2^(1/2)*(s/(x-m))^(1/2)) end proc

 

proc (N) `~`[:-`+`](m, ` $`, `~`[:-`/`](s, ` $`, `~`[:-`^`](`~`[Statistics:-Quantile](Normal(0, 1), Statistics:-Sample(Uniform(1/2, 1), N)), ` $`, 2))) end proc

(4)

Levy := Distribution(
                      PDF=(x->pdf(x)),
                      CDF= (x->cdf(x)),
                      RandomSample = (N->samplingMethod(N))
                    ):

X := RandomVariable(Levy):
PDF(X, x);
print();

# here is an inelegant workaround

ms := {m=0, s=1};
subs(ms, PDF(X, x));
 

piecewise(x < 0, 0, (1/2)*sqrt(2)*sqrt(s/Pi)*exp(-(1/2)*s/(x-m))/(x-m)^(3/2))

 

NULL

 

ms := {m = 0, s = 1}

 

piecewise(x < 0, 0, (1/2)*2^(1/2)*(1/Pi)^(1/2)*exp(-(1/2)/x)/x^(3/2))

(5)

u := Sample(X, 2);

# The substitution doesn't work

subs(ms, u);

u := Vector[row](2, {(1) = m+HFloat(1.7545213270590625)*s, (2) = m+HFloat(0.332342758403606)*s})

 

Vector[row]([m+1.75452132705906*s, m+.332342758403606*s])

(6)

# just a copy-paste of lprint(u): the substitution does work


lprint(u);

subs(ms, Vector[row](2, {1 = m+HFloat(2.49542958430755446)*s, 2 = m+HFloat(4.80772083491271562)*s}, datatype = anything, storage = rectangular, order = Fortran_order, shape = []))

Vector[row](2, {1 = m+HFloat(.529024707466874800)*s, 2 = m+HFloat(1.50334452451578438)*s}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])

 

Vector[row]([2.49542958430755, 4.80772083491272])

(7)

 


 

Download My-own-random-variable.mw

Hi, 

Could you explain me where the third output comes from?

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

Thanks in advance

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

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

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

3 4 5 6 7 8 9 Page 5 of 11