Advanced Search

Question: How to build your own random variable?

Question:How to build your own random variable?

mmcdara 6149 Maple
randomvariable statistics distribution
March 10 2019
1

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
Please Wait...