mmcdara

640 Reputation

10 Badges

3 years, 43 days

MaplePrimes Activity


These are answers submitted by mmcdara

No, Maple hasn't

What are you looking for?

If you just need sampling a Levy(m, s)  distribution (m = location parameter, s = dispersion parameter [classical stat/prob notations], see for instance https://en.wikipedia.org/wiki/Lévy_distribution where m is noted mu and s is noted c) you can use this:

(a few extra material is given to help you define a LevyDistribution)


 

restart:

with(Statistics):

U := RandomVariable(Uniform(1/2, 1)):

N  := 10^4:
SU := Sample(U, N);

SU := Vector(4, {(1) = ` 1 .. 10000 `*Vector[row], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(1)

m := 3;
s := 2;

3

 

2

(2)

SL := m +~ s /~ ( Quantile~(Normal(0, 1), SU) )^~2 ;

SL := Vector(4, {(1) = ` 1 .. 10000 `*Vector[row], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(3)

Histogram(log[10]~(SL), frequencyscale=absolute, minbins=100)

 

For comparison, here is the result returned by R with N=10^4, m=3 and s=2

#

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);

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

(4)

Levy := Distribution(PDF=(x->pdf(x)), CDF= (x->cdf(x)), Conditions = [s > 0])

_m4469968896

(5)

X := RandomVariable(Levy)

_R11002

(6)

Levy__mean := Mean(X);

int((1/2)*_t1*2^(1/2)*(s/Pi)^(1/2)*exp(-(1/2)*s/(_t1-m))/(_t1-m)^(3/2), _t1 = 0 .. infinity)

(7)

Levy__median := Median(X);

(1/2)*(2*RootOf(2*erfc(_Z)-1)^2*m+s)/RootOf(2*erfc(_Z)-1)^2

(8)

evalf(subs({m=0, s=1}, Levy__median));

# compare with wikipedia

evalf(subs({m=0, s=1}, s/2*(1/erfc(1/2))^2));

2.198109338

 

2.174665977

(9)

# plot of the pdf (wiki ways the mode, in case m=0 equals s/3 =1/3 here)

plots:-display(
  plot(subs({m=0, s=1}, PDF(X, x)), x=0..2, gridlines=true),
  plot([[1/3, 0], [1/3, 0.5]], color=blue)
);

 

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

MyLevy := Distribution(PDF=(x->subs(ms, pdf(x))), CDF= (x->subs(ms, cdf(x))))

_m4469962720

(10)

MyX := RandomVariable(MyLevy):

A := Sample(MyX, 10^4, method = [envelope, range = 0.00001..6]):

plots:-display(
   Histogram(A, minbins=100, color=blue, transparency=0.5, view=[0..2, default]),
   plot(subs({m=0, s=1}, PDF(X, x)), x=0..2, color=red, thickness=3,
        title="A scaling problem I've not solved yet\n", font=[Roman, bold, 14])
)

 

 


 

Download Sampling-Levy_Distribution.mw

Here are a few examples.


PS: All RandomVariables have floating point outcomes. If you want to generate "true integer" outcomes, use random instead ... but you will loss the concept of random variable.

RollingDice.mw
 

restart:

with(Statistics):

Dice1 := RandomVariable(DiscreteUniform(1, 6));
Dice2 := RandomVariable(DiscreteUniform(1, 6));

_R

 

_R0

(1)

# example 1:

f1 := (Dice1, Dice2) -> Dice1 + Dice2:
Sample(f1(Dice1, Dice2), 10^5):
Histogram(%, title=f1('Dice1', 'Dice2'));

 

# example 2:

f2 := (Dice1, Dice2) -> Dice1^2 + Dice2:
Sample(f2(Dice1, Dice2), 10^5):
Histogram(%, title=f2('Dice1', 'Dice2'));

 

# example 3:

f3 := (Dice1, Dice2) -> exp(Dice1)*cos(Dice2):
Sample(f3(Dice1, Dice2), 10^5):
Histogram(%, title=f3('Dice1', 'Dice2'));

 

# Finally: a lot of symbolic computations can be done, for instance

print~([seq(Mean__f||n, n=1..3)] =~ Mean~([seq(f||n(Dice1, Dice2), n=1..3)]))
 

Mean__f1 = 7

 

Mean__f2 = 56/3

 

Mean__f3 = (1/36)*exp(1)*cos(1)+(1/36)*exp(2)*cos(1)+(1/36)*exp(3)*cos(1)+(1/36)*exp(4)*cos(1)+(1/36)*exp(5)*cos(1)+(1/36)*exp(6)*cos(1)+(1/36)*exp(1)*cos(2)+(1/36)*exp(2)*cos(2)+(1/36)*exp(3)*cos(2)+(1/36)*exp(4)*cos(2)+(1/36)*exp(5)*cos(2)+(1/36)*exp(6)*cos(2)+(1/36)*exp(1)*cos(3)+(1/36)*exp(2)*cos(3)+(1/36)*exp(3)*cos(3)+(1/36)*exp(4)*cos(3)+(1/36)*exp(5)*cos(3)+(1/36)*exp(6)*cos(3)+(1/36)*exp(1)*cos(4)+(1/36)*exp(2)*cos(4)+(1/36)*exp(3)*cos(4)+(1/36)*exp(4)*cos(4)+(1/36)*exp(5)*cos(4)+(1/36)*exp(6)*cos(4)+(1/36)*exp(1)*cos(5)+(1/36)*exp(2)*cos(5)+(1/36)*exp(3)*cos(5)+(1/36)*exp(4)*cos(5)+(1/36)*exp(5)*cos(5)+(1/36)*exp(6)*cos(5)+(1/36)*exp(1)*cos(6)+(1/36)*exp(2)*cos(6)+(1/36)*exp(3)*cos(6)+(1/36)*exp(4)*cos(6)+(1/36)*exp(5)*cos(6)+(1/36)*exp(6)*cos(6)

 

[]

(2)

 

 

 

Like Tom I'm not sure to correctly understand your problem.
It seems to me you want to sample truncated Poisson distributions (between 0 and 10 included) whose parameter lambda is itself a linear combination of random variables (am I right? )

If yes, this piece of code will possibly answer your question.
Be careful: not all the values of lambda you create from your model are positive!!!
In the attached file I show you how to estimate the frequency of negative lambda values. This frequency goes to zero as the value of N (that is the size of the vector lambda) tends to +infinity. For N=10 the ratio of negative lambda values is close to 1/4.


 

restart:

with(Statistics):

randomize():

N    := 10;
x__0 := Vector[row](P, [1$N]);
x__1 := Sample(Binomial(N, 0.4), N);
x__2 := Sample(Normal(0, 1), N);

N := 10

 

`#msub(mi("x"),mi("0"))` := Vector[row](10, {(1) = 1, (2) = 1, (3) = 1, (4) = 1, (5) = 1, (6) = 1, (7) = 1, (8) = 1, (9) = 1, (10) = 1})

 

`#msub(mi("x"),mi("1"))` := Vector[row](10, {(1) = 5.0, (2) = 1.0, (3) = 6.0, (4) = 3.0, (5) = 5.0, (6) = 2.0, (7) = 7.0, (8) = 2.0, (9) = 6.0, (10) = 8.0}, datatype = float[8])

 

`#msub(mi("x"),mi("2"))` := Vector[row](10, {(1) = .5691573712734501, (2) = .10685006448894822, (3) = -.39789312201697974, (4) = -.5395405825678521, (5) = .6183149648978897, (6) = .5512360714665099, (7) = 1.2991554584043454, (8) = .18083457234126987, (9) = -.47061639469370886, (10) = 0.8967688506232983e-1}, datatype = float[8])

(1)

local lambda:

lambda := unapply(add(beta__||k *~ x__||k, k=0..2), (seq(beta__||k, k=0..2)) );

 

 

B := (alpha, NMC) -> Sample(Poisson(alpha), NMC, method=[discrete, range=0..10])

proc (alpha, NMC) options operator, arrow; Statistics:-Sample(Poisson(alpha), NMC, method = [discrete, range = 0 .. 10]) end proc

(3)

beta := (0.25, 0.5, 3.2);
lambda~(beta);

beta := .25, .5, 3.2

 

Vector[row]([4.57130358807504, 1.09192020636463, 1.97674200954566, 0.234701357828733e-1, 4.72860788767325, 3.01395542869283, 7.90729746689391, 1.82867063149206, 1.74402753698013, 4.53696603219946])

(4)

MC := 10:

NumberOfNegativeLambdaValues := 0:

U := NULL:

for u in lambda~(beta) do
   if u > 0 then
     b := round~(convert(B(u, MC), list)) ;
     U := U, b[];
   else
     # printf("Error, parameter is negative(%a)\n", u):
     NumberOfNegativeLambdaValues := NumberOfNegativeLambdaValues + 1:
   end if;
end do;

printf("%2.f %% of the lambda have a negative value\n", evalf(NumberOfNegativeLambdaValues/MC)*100);

DataSummary([U]);
Histogram([U], discrete, minbins=N+1, thickness=20)

 0 % of the lambda have a negative value

 

Vector[column]([[mean = 3.30000000000000], [standarddeviation = 2.75790873780490], [skewness = .542012393898590], [kurtosis = 2.20622450437294], [minimum = 0.], [maximum = 10.], [cumulativeweight = 100.]])

 

 

REMARKS

The frequency of occurrences of negative lambda values can be easily determined

RV__0 := 1:
RV__1 := RandomVariable(Binomial(N, 4/10));
RV__2 := RandomVariable(Normal(0, 1));

RV__3 := unapply( add(beta__||k *~ RV__||k, k=0..2), (seq(beta__||k, k=0..2)) )

_R11

 

_R12

 

proc (beta__0, beta__1, beta__2) options operator, arrow; _R11*beta__1+_R12*beta__2+beta__0 end proc

(5)

(Mean, StandardDeviation)(RV__3(seq(beta__||k, k=0..2)) );

4*beta__1+beta__0, (1/5)*(60*beta__1^2+25*beta__2^2)^(1/2)

(6)

(Mean, StandardDeviation)(RV__3(beta) );

FrequencyOfNegativeLambda := Probability(RV__3(beta) <= 0);

2.25, 3.292415527

 

.2472407931

(7)

 


 

Download Poisson.mw


 

Statistics:-LinearFit can do the job quite well. So there is no bug in Statistics:-LinearFit  (or at least not on this problem).

The issue comes from the fact that the problem is not relevant, per se, of the framework of Ordinary Least Squares (OLS) regression where we aim to assess the values of the parameters A and B of the model:

" Conditional expectation of Y given X=x = B*x+A "

The reason for this is that the slope is known.
In the attached file I show how a rewritting of the initial problem sets it back in the OLS framework. But this is neither immediate not natural.

However, the simplest method to assess the value of "a" comes from classical estimation results in OLS : the best unbiased estimator of "a" is just the mean of pts2 - 3*pts1

 


 

restart:

with(Statistics):

pts1 := Vector([3, 5, 21]);

pts2 := Vector([4, 6, 16])

pts1 := Vector(3, {(1) = 3, (2) = 5, (3) = 21})

 

pts2 := Vector(3, {(1) = 4, (2) = 6, (3) = 16})

(1)

It couldn't be easier

# From standard results in Ordinary Least Squares (OLS) regression,
# the Best Linear Unbiased Estimator (BLUE) of the intercept A
# in the model: Expectation(Y | X) = B*X+A is just
# BLUE(A) = Mean(Y) - BLUE(B)*Mean(X)
#
# Here BLUE(B) is supposed to be known and set to 3.
# Then

BLUE__A := Mean(pts2 - 3 *~ pts1);
add (pts2 - 3 *~ pts1) / 3

HFloat(-20.333333333333332)

 

-61/3

(2)

Using Statistics:-LinearFit

# The problem is not relevant of OLS but can be solved with LinearFit
# (at the price of a few pirouettes)
#
# Rewrite the problem this way:
#   1/ define Z as Y-3*X

Z := (X, Y) -> Y - 3*X;

#   2/ imagine you repeat 3 times the same experiment for a same value of the
#      control variable "a" ("a" is the only regressor of the problem)
#      Write Newpts2 the vector of the outcomes of these experiments
Newpts2 := Z~(pts1, pts2);
#   3/ define then the "design matrix" as an arbitrary constant vector,
#      for instance < 1, 1, 1>

Newpts1 := Vector([1, 1, 1]);

#   4/ apply LinearFit as usual, remembering that the only regressor is "a"

L := LinearFit([a], Newpts1, Newpts2, a);

#   5/ Substitute "a" by the value it takes in the vector Newpts1

A := subs(a=1, L);
convert(%, rational);

Z := proc (X, Y) options operator, arrow; Y-3*X end proc

 

Newpts2 := Vector(3, {(1) = -5, (2) = -9, (3) = -47})

 

Newpts1 := Vector(3, {(1) = 1, (2) = 1, (3) = 1})

 

-HFloat(20.33333333333334)*a

 

HFloat(-20.33333333333334)

 

-61/3

(3)

 

 


 

Download Simplest.mw

Simulation of a "true" Galton board is very difficult.
You have, in all rigor, to simulate the drop of a sphere (a disc in 2D) shocking many fixed cylindrical pins (discs in 2D).
The trajectory of a disc (assuming the problem is 2D) is described by a succession of free flights submitted to the gravitational field (rather simple to code), each of them ending when the falling disc hurts a fixed pin (another disc).
Determining the impact point is a difficult problem.
More of this you have to describe the physics of the shock (elastic or not), account for possible slipping and rolling...
Last but not least a true Galton board has a limited width an the pins are confined in a kind of funnel whose walls modify the fall at at boundaries of the board.

What I propose you is a simpler simulation based on a random walk on a "triangular lattice" (Pascal's lattice).
Starting from the position (i, j) a "ball" can move either to position (i+1, j+1) or (i+1, j-1) with equal probabilities. The process is repeated from the initial position (0, 0) a given number of times.
Think to this to some kind of "discrete Galton board".

Look to the animation in the attached file.
If it suits you I will give you some explanation about it... if required


WATCH OUT: I faced difficulties to load the file with the animations, so I commented the corresponding line (red written on yellow background).
To run the animation, once the plot after the command Fallouts(100, 10) has been displayed, just clik on it and run the animation by using the keys in the upper left of your Maple worksheet.

It's likely that this code can be written in clever and more performing way.
 

restart:

see for instance https://en.wikipedia.org/wiki/Random_walk

with(Statistics):
with(plots):
with(plottools):

randomize():

R := rand(0 .. 1):

Fallouts := proc(NP, N)
   local  Fallout:
   local  GaltonBoard, depart, t, n, c, T, GBH;
   global cs, GB:
   
 
   # fallout of a single "ball"
 

   Fallout := proc(np)
      local GBH, depart, t, n, c, T;
      global cs, GB:


      GBH := copy(GB):
      depart := [0, 0];
      t      := depart:
      for n from 1 to N do
         depart := depart +~ [1, (2*R()-1)];
         t      := t, depart
      end do:   
      t  := [t, [depart[1]+1, depart[2]]];
      c  := unapply(CurveFitting:-Spline(op~(1, t), op~(2, t), x, degree=1), x);
      T  := transform((x,y) -> [y+N+1, x]):
  
      if np > 1 then
         GBH := display(GBH, T(Histogram([cs], discrete, thickness=5, frequencyscale=absolute)) )
      end if;
      
      cs := cs, depart[2];
   
   plots[animate](plot, [c(x), x=0..A, color=red, thickness=2], A=0..N+1, frames=N+2, background=GBH);
   
   end proc:

   GaltonBoard := [seq(seq([i, j], j in [seq(-i..i, 2)]), i=0..N)]:
   GB := plot(GaltonBoard, style=point, symbol=solidcircle, symbolsize=10, scaling=constrained, gridlines=true, color=blue):

   cs := NULL:

   plots[animate]( Fallout, [np], np=1..NP, frames=NP);

end proc:

global cs:

# drop 100 "balls" on a "discrete Galton Board" made oh 10 layers of pins

# Fallouts(100, 10);   # TO BE UNCOMMENTED

cs

cs

(1)

# Histogram([cs])

 


 

Download Galton-RandomWalk.mw

Maybe you could try this


 

restart

with(Statistics):

XG := h -> RandomVariable(Normal(0, h))

proc (h) options operator, arrow; Statistics:-RandomVariable(Normal(0, h)) end proc

(1)

MaxIt := 1000

1000

(2)

randomize():

# Be careful: h must be strictly positive

h := evalf([seq(1-k/MaxIt, k=0..MaxIt-1)]):
h[1],h[-1];

1., 0.1000000000e-2

(3)

# do you really need to print the results ?
#
# for i from 1 to MaxIt do
#    printf("%a\n%a\n", Sample(XG(h[i]), 5), i);
# end do:

# Maybe this could suit you?


# Stage 1 ; generate a random matrix by sampling Normal(0, 1)

M := Sample(XG(1), [MaxIt, 5], method = [envelope, range = -1 .. 1]):


# Watch out: stage 2 is computationnaly expensive.
#            It is presented here to help understanding the method
#
# Stage 2: rescale M according to the standard deviations 1, 0.999, ...0.001

S := LinearAlgebra:-DiagonalMatrix(h):

# Stage 3: do the product S.M

Res := S.M

Res := Vector(4, {(1) = ` 1000 x 5 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(4)

# The same thing in a more efficient way

M   := CodeTools:-Usage( Sample(XG(1), [MaxIt, 5], method = [envelope, range = -1 .. 1]) ):
Res := CodeTools:-Usage( <  seq(M[n, ..] *~ h[n], n=1..MaxIt) > );

memory used=0.73MiB, alloc change=0 bytes, cpu time=35.00ms, real time=11.00ms, gc time=0ns
memory used=0.60MiB, alloc change=0 bytes, cpu time=9.00ms, real time=3.00ms, gc time=0ns

 

Res := Vector(4, {(1) = ` 1000 x 5 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(5)

for i from 1 to MaxIt do
  # printf("%3d : %a\n", i, [entries(Res[i,..], nolist)]);
end do:

 


 

Download Gaussian_Reply.mw

Several people have already answered your question in different ways

To complete their propositions, here is a time and memory comparison of different ways to answer yout problem

restart:

with(Statistics):   # the most natural (?) and one of the fastest if you want 0 and 1 outputs
B := RandomVariable(Bernoulli(1/2)):

numThrows:= 10^5:

CodeTools:-Usage( Sample(B, numThrows) ):

memory used=0.87MiB, alloc change=0 bytes, cpu time=7.00ms, real time=7.00ms, gc time=0ns

 

restart;

with(StringTools) : # a simple way to get "H" and "T" ouputs

numThrows:= 10^5:
CodeTools:-Usage( Random(numThrows, "HT" ) ):                     # a string of the form "HHT..."
CodeTools:-Usage( Explode(Random(numThrows, "HT" )) ):            # a list of the form [ "H", "H", "T", ...
CodeTools:-Usage( Vector(Explode(Random(numThrows, "HT" ))) ):    # a vector (here fot comparison with tomleslie's solution below)

memory used=154.93KiB, alloc change=0 bytes, cpu time=193.00ms, real time=193.00ms, gc time=0ns
memory used=0.86MiB, alloc change=0 bytes, cpu time=192.00ms, real time=193.00ms, gc time=0ns

memory used=22.27MiB, alloc change=0 bytes, cpu time=383.00ms, real time=384.00ms, gc time=41.38ms

 

restart;

# tomleslie  ... probably the fastest (?) to deliver "H" and "T" ouputs
numThrows:= 10^5:

r:= rand(1..2):

CodeTools:-Usage(  Vector[row](numThrows, i->r()) ):
CodeTools:-Usage(  Vector[row](numThrows, i->["H","T"][r()]) ):

memory used=0.81MiB, alloc change=0 bytes, cpu time=33.00ms, real time=33.00ms, gc time=0ns
memory used=0.77MiB, alloc change=0 bytes, cpu time=38.00ms, real time=38.00ms, gc time=0ns

 

restart;

with(LinearAlgebra): # advantage: there exist a very fast way to tansform 0 & 1 into "H" & "T"

numThrows:= 10^5:

V := CodeTools:-Usage(  RandomVector[row](numThrows,generator=rand(0..1)) ):  # outputs are 0 or 1
W := CodeTools:-Usage(  subs({0="H", 1="T"}, V) ):                            # outputs are "H" or "T"

memory used=4.70MiB, alloc change=2.00MiB, cpu time=50.00ms, real time=51.00ms, gc time=5.63ms
memory used=0.76MiB, alloc change=0 bytes, cpu time=3.00ms, real time=3.00ms, gc time=0ns

 

restart;

with(SignalProcessing):  # comparable to Statistics

numThrows:= 10^5:

V := CodeTools:-Usage( floor~(GenerateUniform(numThrows, 0, 2)) ):  # the fastest (?), use Hfloats (default)

UseHardwareFloats := false:
CodeTools:-Usage( floor~(GenerateUniform(numThrows, 0, 2)) ):  # don't use Hfloats

memory used=1.79MiB, alloc change=0 bytes, cpu time=7.00ms, real time=8.00ms, gc time=0ns

memory used=203.74MiB, alloc change=47.02MiB, cpu time=1.58s, real time=1.59s, gc time=88.87ms

 

restart:

with(combinat):  # just for fun

numThrows:= 10^5:

CodeTools:-Usage( randperm(["T"$numThrows, "H"$numThrows])[1..numThrows] ):

memory used=20.67MiB, alloc change=4.58MiB, cpu time=161.00ms, real time=162.00ms, gc time=61.65ms

 

 


 

Download Still_Another_Way.mw

 

After some digressions unrelated to your initial question, here is my definitive contribution

final.mw

restart:

with(Statistics):

Seed := randomize():

alias(RV=RandomVariable):

# Number of dice

N := 3;

3

(1)

# These two lines can be replaced by any algorithm amid those presented in this
# same thread (just choose the one you want)


S := add(RV(DiscreteUniform(1, 6)), k=1..N):
P := Probability~(S =~ [$N..6*N])

[1/216, 1/72, 1/36, 5/108, 5/72, 7/72, 25/216, 1/8, 1/8, 25/216, 7/72, 5/72, 5/108, 1/36, 1/72, 1/216]

(2)

# Define a discrete RV with outcomes of probabilities given by P
#
# Watch out: by default (I do not know if and how this could be changed), this RV has outcomes 1..5*N+1
# This raises an issue that will be fixed when the resuklts are ploted (the outcomes are N..6*N)

Dice := RV(ProbabilityTable(P))

_R2

(3)

# Sample "Dice" as any other RV

M := 100:
X := Sample(Dice, M):

h := Histogram(X, minbins=5*N+1):

# The next operation is aimed to superimpose the exact mss function (aka P) and its
# empirical estimation h.
#
# I transform h into a list that can be plot with ColumnGraph.
# Watch out:
#      Maple 2018: use op~(2, op~(1, [op(1..-2, op(1, h))] ));
#      Maple 2015: use map(u -> print(u[-1][2]), [op(1..-2, op(1, h))]);

interface(version);

H := map(u -> u[-1][2], [op(1..-2, op(1, h))]):

ColumnGraph(
    [H, P],
    legend=["Sampled", "Exact"],
    color=[red, yellow],
    background=gray,
    title=cat("Sample size = ", M),
    tickmarks=[[seq(i=i+N, i=0..5*N)], default],
    size=[600, 400]
)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

 

 

 


 

Download final.mw

 

 

@tomleslie 

A few remarks:

  1. Use Quantile(X, k+i+j, numeric) for a faster generation:
    CodeTools:-Usage(Matrix(100, 10, (i,j) -> Quantile(X, (i-1)/100+(j-1)/1000, numeric))):
    memory used=14.25MiB, alloc change=0 bytes, cpu time=158.00ms, real time=159.00ms, gc time=0ns
    CodeTools:-Usage(Matrix(100, 10, (i,j) -> Quantile(X, (i-1)/100+(j-1)/1000))):
    memory used=430.05MiB, alloc change=0 bytes, cpu time=6.28s, real time=5.90s, gc time=555.98ms

    Note the CPU time (in this particular case) can be halved: Quantile(X, p) = -Quantile(X, 1-p)
     
  2. Avoid computing Quantile(X, 0) which is obviously -infinity and not -15.681...
     

Here is a version which uses DocumentTools for a smarter rendering

QuantileTable.mw


 

As you wrote it, your maplet can't work: Evaluate needs a function as argument, not a sequence of instructions.
Here is the correct way to code it in a simple case (no "for a in .... end do" loop ... but you will be able to generalize with from Carl's answer).

PS: the code I wrote conforms to your's, but it could be improved (remark first that there is no need to click on the "Display function" function button to visualize the MathML code in the MMLV field). It's also better to define "Input function  f(x)=" within a Label component, ...

restart:

with(Maplets[Elements]):
with(Student[Calculus1]):

f := proc()
local g, TL2:

g := proc(func, slp)
  local a;
  a := Tangent(func,x=slp):
  Maplets:-Tools:-Set('MMLV1'=a)
end proc:


TL2:=Maplet([
  [
    "Input function  f(x)=", TextField['func'](20)
  ],
  [
    Button("Display function", Action(Evaluate('MMLV'='MathML[Export](func)')))
  ],
  MathMLViewer['MMLV'](),
  [
    "Input slope", TextField['slp'](10)
  ],  
  [
    Button("Tangent line", Action( Evaluate('function'=g, Argument(func), Argument(slp))))
  ],   
  MathMLViewer['MMLV1']()   
]):

Maplets[Display](TL2):

end proc:

f()

 


 

Download Maplet.mw

 

You had forgotten that the projected dx*dy surface wasn't the surface itself:


 

restart;

with(VectorCalculus):

SetCoordinates(cartesian[x, y, z]):

R := sqrt(2):
s := x^2+y^2+z^2-R^2;
u := VectorField([x, 1, z]);

N := Gradient(s);

x^2+y^2+z^2-2

 

u := Vector(3, {(1) = x, (2) = 1, (3) = z}, attributes = [vectorfield, coords = cartesian[x, y, z]])

 

N := Vector(3, {(1) = 2*x, (2) = 2*y, (3) = 2*z}, attributes = [vectorfield, coords = cartesian[x, y, z]])

(1)

n := N/sqrt(add(N[k]^2, k = 1 .. 3));
u . n;

un := algsubs(x^2+y^2+z^2=R^2, %):
un := simplify(un)

n := x*`#mover(mi("e"),mo("&lowbar;"))`[x]/sqrt(x^2+y^2+z^2)+y*`#mover(mi("e"),mo("&lowbar;"))`[y]/sqrt(x^2+y^2+z^2)+z*`#mover(mi("e"),mo("&lowbar;"))`[z]/sqrt(x^2+y^2+z^2)

 

x^2/(x^2+y^2+z^2)^(1/2)+y/(x^2+y^2+z^2)^(1/2)+z^2/(x^2+y^2+z^2)^(1/2)

 

-(1/2)*2^(1/2)*y^2+(1/2)*2^(1/2)*y+2^(1/2)

(2)

intr1 := solve(subs(z = 0, s), y);

intr2 := solve(subs(z = 0, y = 0, s), x);

(-x^2+2)^(1/2), -(-x^2+2)^(1/2)

 

2^(1/2), -2^(1/2)

(3)

# forgot that the projected dx*dy surface is not the surface irself

l     := sqrt(x^2+y^2);
h     := sqrt(2-x^2-y^2);
phi   := arctan(h/l);
scale := simplify(1/sin(phi)) assuming x^2+y^2 > 0;

# surface of the sphere

R1 := 2*int(int(scale, y = intr1[2] .. intr1[1]), x = intr2[2] .. intr2[1], numeric);

evalf(4*Pi*R^2);

(x^2+y^2)^(1/2)

 

(-x^2-y^2+2)^(1/2)

 

arctan((-x^2-y^2+2)^(1/2)/(x^2+y^2)^(1/2))

 

2^(1/2)/(-x^2-y^2+2)^(1/2)

 

25.13274122

 

25.13274123

(4)

R1 := 2*Student:-MultivariateCalculus:-MultiInt(scale*un, y = intr1[2] .. intr1[1], x = intr2[2] .. intr2[1])

(16/3)*2^(1/2)*Pi

(5)

evalf(%);

23.69537567

(6)

du := Divergence(u)

2

(7)

intr1 := solve(s, z);

intr2 := solve(subs(z = 0, s), y);

intr3 := solve(subs(z = 0, y = 0, s));

R2 := int(int(int(du, z = intr1[2] .. intr1[1]), y = intr2[2] .. intr2[1]), x = intr3[2] .. intr3[1]);

evalf(%)

(-x^2-y^2+2)^(1/2), -(-x^2-y^2+2)^(1/2)

 

(-x^2+2)^(1/2), -(-x^2+2)^(1/2)

 

2^(1/2), -2^(1/2)

 

(16/3)*2^(1/2)*Pi

 

23.69537567

(8)

 


 

Download GaussXYZ.mw

In your particular case it's far more easier to verify the Gass theorem in spherical coordinates

I'll give a closer look to your code later
 

restart;

with(IntegrationTools):
with(VectorCalculus):
SetCoordinates(cartesian[x, y, z]):

MySphere := x^2+y^2+z^2=2;
u := VectorField([x, 1, z]);

N := Gradient(lhs(MySphere));

x^2+y^2+z^2 = 2

 

u := Vector(3, {(1) = x, (2) = 1, (3) = z}, attributes = [vectorfield, coords = cartesian[x, y, z]])

 

N := Vector(3, {(1) = 2*x, (2) = 2*y, (3) = 2*z}, attributes = [vectorfield, coords = cartesian[x, y, z]])

(1)

n := N/sqrt(add(N[k]^2, k = 1 .. 3));
n := subs(MySphere, n);

n := x*`#mover(mi("e"),mo("&lowbar;"))`[x]/sqrt(x^2+y^2+z^2)+y*`#mover(mi("e"),mo("&lowbar;"))`[y]/sqrt(x^2+y^2+z^2)+z*`#mover(mi("e"),mo("&lowbar;"))`[z]/sqrt(x^2+y^2+z^2)

 

n := Vector(3, {(1) = (1/2)*2^(1/2)*x, (2) = (1/2)*2^(1/2)*y, (3) = (1/2)*2^(1/2)*z}, attributes = [vectorfield, coords = cartesian[x, y, z]])

(2)

un := DotProduct(u, n);
du := Divergence(u);

(1/2)*x^2*2^(1/2)+(1/2)*2^(1/2)*y+(1/2)*z^2*2^(1/2)

 

2

(3)

SphereRadius := sqrt(rhs(MySphere));

2^(1/2)

(4)

J := (r, phi) -> r^2*sin(phi);

int(int(int(du*J(r, phi), theta=0..2*Pi), phi=0..Pi), r=0..SphereRadius);
evalf(%);

proc (r, phi) options operator, arrow; VectorCalculus:-`*`(r^2, sin(phi)) end proc

 

(16/3)*Pi*2^(1/2)

 

23.69537567

(5)

int(int(subs(r=SphereRadius, unsph)*J(SphereRadius, phi), theta=0..2*Pi), phi=0..Pi);
evalf(%);

(16/3)*Pi*2^(1/2)

(6)

 


 

Download Gauss.mwGauss.mw

Not the solution, just a way to derive the equations by applying differential operators.
I think you have an ambiguity about the pressure (you declare p bar as a constant but it has a priori different values at z=-1 and z=+1 (unless pa is null)

For the solution: I'm not familiar with pdsolve and PDETools, but I'm afraid that you will have to do the things by hand. 
If I'm not mistaken this is a classical poiseuille flow with uniform velocity at the inlet (parabolic profile expected at the outlet) and I remember solving it by hand long ago at school. Maybe you will have to reproduce explicitely the same modus operandi in Maple?

Wait for a more enlightened reply about pdsolve/PDEtools


 

restart:

with(Student[VectorCalculus]):

SetCoordinates(cylindrical[r,theta, z]):

alias(u=u(r, 'theta', z));
alias(v=v(r, 'theta', z));
alias(p=p(r, 'theta', z));

u

 

u, v

 

u, v, p

(1)

V := VectorField(< v, 0, u >);
Divergence(V);

V := Vector(3, {(1) = v(r, theta, z), (2) = 0, (3) = u(r, theta, z)}, attributes = [vectorfield, coords = cylindrical[r, theta, z]])

 

(v+r*(diff(v, r))+r*(diff(u, z)))/r

(2)

Continuity_equation := eval(expand(%));

v/r+diff(v, r)+diff(u, z)

(3)

Laplacian(V):
Diffusion_term := beta^2 *~ map(expand, eval(%));

Diffusion_term := Vector(3, {(1) = beta^2*(-v(r, theta, z)/r^2+(diff(v(r, theta, z), r))/r+diff(diff(v(r, theta, z), r), r)+(diff(diff(v(r, theta, z), theta), theta))/r^2+diff(diff(v(r, theta, z), z), z)), (2) = 2*beta^2*(diff(v(r, theta, z), theta))/r^2, (3) = beta^2*(diff(diff(u(r, theta, z), z), z)+(diff(u(r, theta, z), r))/r+diff(diff(u(r, theta, z), r), r)+(diff(diff(u(r, theta, z), theta), theta))/r^2)})

(4)

Pressure_term := 1 *~ Gradient(p)

Pressure_term := Vector(3, {(1) = diff(p(r, theta, z), r), (2) = (diff(p(r, theta, z), theta))/r, (3) = diff(p(r, theta, z), z)})

(5)

Advection_term := 1 *~ Vector[column](3, i -> DotProduct(V, Gradient~(V[i])))

Advection_term := Vector(3, {(1) = v(r, theta, z)*(diff(v(r, theta, z), r))+u(r, theta, z)*(diff(v(r, theta, z), z)), (2) = 0, (3) = v(r, theta, z)*(diff(u(r, theta, z), r))+u(r, theta, z)*(diff(u(r, theta, z), z))})

(6)

C := :-Vector[column](3, i -> Advection_term[i] = - Pressure_term[i] + Diffusion_term[i])[[1, 3]]

C := Vector(2, {(1) = v(r, theta, z)*(diff(v(r, theta, z), r))+u(r, theta, z)*(diff(v(r, theta, z), z)) = -(diff(p(r, theta, z), r))+beta^2*(-v(r, theta, z)/r^2+(diff(v(r, theta, z), r))/r+diff(diff(v(r, theta, z), r), r)+(diff(diff(v(r, theta, z), theta), theta))/r^2+diff(diff(v(r, theta, z), z), z)), (2) = v(r, theta, z)*(diff(u(r, theta, z), r))+u(r, theta, z)*(diff(u(r, theta, z), z)) = -(diff(p(r, theta, z), z))+beta^2*(diff(diff(u(r, theta, z), z), z)+(diff(u(r, theta, z), r))/r+diff(diff(u(r, theta, z), r), r)+(diff(diff(u(r, theta, z), theta), theta))/r^2)})

(7)

NS := [Continuity_equation, entries(C, nolist)]

[v/r+diff(v, r)+diff(u, z), v*(diff(v, r))+u*(diff(v, z)) = -(diff(p, r))+beta^2*(-v/r^2+(diff(v, r))/r+diff(diff(v, r), r)+(diff(diff(v, theta), theta))/r^2+diff(diff(v, z), z)), v*(diff(u, r))+u*(diff(u, z)) = -(diff(p, z))+beta^2*(diff(diff(u, z), z)+(diff(u, r))/r+diff(diff(u, r), r)+(diff(diff(u, theta), theta))/r^2)]

(8)

# You can say "p bar is a constant and { p(z=-1)=0 & p(z=-1)=0 & p(z=-+1)=p__a } unless p__a=0 too  

BC := [
:-D[1](u)(0, 'theta', z) = 0,
v(0, 'theta', z) = 0,
u(1, 'theta', z) = 0,
v(1, 'theta', z) = p+alpha
]:

seq(print(BC[i]), i=1..numelems(BC))

(D[1](u))(0, theta, z) = 0

 

v(0, theta, z) = 0

 

u(1, theta, z) = 0

 

v(1, theta, z) = p+alpha

(9)

 


 

Download NavierStokes.mw

Replace your line   Diam;
by                           Diam := simplify(Diam);

Nevertheless I would be tempted to say that Maple should simplify automatically.

Your line evalf(Diam) is corrupted.
Just rewrite it

TelescopeUnits_reply.mw


PS:Personnaly I never use this document mode input that I find very capricious and subjected to hide characters mistakenly typed

I'm not sure I understand you perfectly.
What I think is: you have a discrete RV "X" with mass function fX(x)=1/(x+1) for x = 0, ..., k, where k is some given number (>0).

If It is that, a way to proceed is to use a ProbabilityTable Distribution

 

restart:

with(Statistics):

# Discrete distribution from a probability table
# The sum of all masses must be equal to 1.

P := k -> [seq(1/(x+1), x=0..k)] /~ add(1/(x+1), x=0..k) ;

# example:

K := 5:
N := 10^3:

P(K);

X := RandomVariable(ProbabilityTable(P(K))):

S := Sample(X, N):
S[1..10];
Histogram(S);

 

P := proc (k) options operator, arrow; `~`[:-`/`]([seq(1/(x+1), x = 0 .. k)], ` $`, add(1/(x+1), x = 0 .. k)) end proc

 

"[[Typesetting:-mfrac(Typesetting:-mn("20",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("10",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("20",mathvariant = "normal"),Typesetting:-mn("147",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("5",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("4",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("10",mathvariant = "normal"),Typesetting:-mn("147",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false")]]"

 

Vector[row]([3., 2., 6., 2., 1., 1., 2., 1., 2., 1.])

 

 

 


 

Download DiscreteRV.mw

2 3 4 5 6 7 8 Page 4 of 9