535 Reputation

11 Badges

8 years, 189 days

MaplePrimes Activity

These are replies submitted by sand15


It's strange you get Matrix(3, 3, {(1, 1) = -1, (1, 2) = -1, (1, 3) = -1, (2, 1) = -2, (2, 2) = -2, (2, 3) = -2, (3, 1) = -3, (3, 2) = -3, (3, 3) = -3}) with the command Matrix(3, 3, `-`);

In my case (Maple2018, worksheet mode) I get

Matrix(%id = 18446746664292233078), just as with Matrix(3, 3, (i,j) -> i-j);


The result returned by  Matrix(3, 3, `.`); is very curious, would you have any idea about it?

@Carl Love 

Right Carl, good point (beta distribution) and nice work (the completely formal solution).

It would be interesting to try to generalize this formal derivation to non uniform pdf.

For some of them order distribution are known and I think it could be possible to obtain interesting results.

By the way, my last mistake was that I write Ecost as a sum of two terms invoking the integral of Q-x as I should have written that it was defined as the value work Wich these two integrals are equal.

I think it would be better to avoid integrating Q-x between a and Q and x-Q between Q and b: fundamentally this is the same loss function Q-x (or x-Q) that we use.


Hope it's clear, I send this from my smartphone


@cinderella look to my code and your code.

The only difference is the expression of N. I define N as the integral of Q-x against the measure f, as you define  this same N as the integral of x-Q against the measure f.

If you are sure of what you did change the sign of N in my code.

You can do this in the 'numerical' code l sent you yesterday or in the formal code I sent you nine hours ago (in this matter change the sign of MAX in the expression of Qstar) 

I have no doubt you will be capable to fix this by yourself. Now I'm in weekend, see you in two days if necessary.


Then you can use this
reference : Maple V, release 5, Monagan et al, 1996  page 55

uniform := proc(r::constant..constant)
  local intrange, f:
  intrange := map(x -> round(x*10^Digits), evalf(r) ):
  f := rand(intrange)
  (evalf @ eval(f)) / 10^Digits)
end proc:

# usage

N := 100:
U := uniform(1..N):
s := [seq(U(), n=1..N)]:
t := sort(s, output=permutation)


This is indeed a fine way to address the problem.
Personally, despite an intensive use of the Statistics package, I'd never payed attention to this EmpiricalDistribution feature.
Maybe you have given a definitive answer to the question, I'm interested in reading acer's opinion on this point.




(mmcdara from my workplace)

"that integration is (for this kind of example) an alternative" : right, is an adhoc alternative limited to this special case of a function of 2 dice.

"But the full set of inputs -- to be passed to g -- is still produced with the nested seq. (The RVs assigned to Dice1 and Dice2 are not used.)"  : I assume you're referring to this instruction?

outcomes := [ { seq(seq(g(u,v), u in [$1..6]), v in [$1..6]) }[] ];

Right again, the whole construction I did is a liittle bit artificial and very specific to the problem.

More if this I think it would be better to be able to compute the probability for every integer values between -4 and 4 included.


Je ne suis pas certain que demander à ce qu'on vous fasse vos TD/TP de licence ou d'école d'ingé soit une bonne manière de vous faire apprendre Maple ?
Vous auriez pu, au moins ne pas reproduire mot à mot le sujet de ce TP.

I'm not sure that asking for someone to do your master's or  engineering school's tutorials (or practical work) is a good way to get you to learn Maple?
You could have at least not reproduced word for word the subject of your practical work..

As a starting point, for question Aa:
First of all you need a stopping criterion, for instance a measure of the difference between the exact value (if known) and its iterated approximations x[n], or the difference between x[n] and x[n-1].
Here is an example :

babylon := proc(x0, k, epsilon)
  local x:
  x := copy(x0):
  while abs(x^2-k) > epsilon do
    x := (x+k/x)/2
  end do:
  return x;
end proc:

# usage
babylon(2, 1, 0.001)


For point Ab: search list, or vector, or Array, ... in the help pages

For point B: search "thisproc" in the help page and take inspiration from the recursive calculus of the factorial of an integer


Of course!

I'm stupid, I didn't even think to use Sample~ !!!

Thanks Tom


Just read the help pages about Grid:-Launch.
You will find an example (primetest) that  I've just transposed to your test case.

Unfortunately I can give you detailed explanations, beyond that of saying "respect the content of the help pages".
If you are interesting in knowing why you must write babyGrid := proc() instead of babyGrid := proc(X) and later Grid:-Launch(babyGrid, imports=['X']) instead of Grid:-Launch(babyGrid, X), I'm not qualified to answer.

All you have to do is wait for someone else to give you this information.

@Carl Love  @tomleslie


Before sending my answer to Erik, I had tried to sample a function of two RV with values in R^2 (typically the "f" function used by Tom).

Remember I used RandomVariables Dice 1 and Dice2.
So I wrote, in a natural way,  f := (Dice1, Dice2) -> [Dice1, Dice2] and Sample(f(Dice1, Dice2), N).
As it generated an error, I tried to lure Maple by writting instead  f := (Dice1, Dice2) -> Dice1+I*Dice2
... which generated the same error as before.

All of this seems related to the impossibility (?) to define a multivariate random variable.
Am I right ?
Is it possible (without sampling independantly Dice1 and Dice 2 as Tom did) to define multivariate RV in Maple?



For as far as  I know the only limitation to Excel (pack office 12) is that the number of rows must not exceed 32767 (maybe 32768?)
I guess the same kind of limitations applies on the number of columns too, this regardless to the Excel version.
As you evoked R, my advice is to use ExportMatrix instead of ExcelTools:-Export or your manual procedure (it's always better to do the job programmatically when possible)

@acer @tomleslie

What version are you using?

  • It was Maple 2016 (usually I use Maple 2018 but I had have a long simulation running on with it and, in these case, I prefer to open a new worksheet in an older version).
    I guess that If I had done the test with Maple 2018 I would not sent this question

Does it change if you set interface(typesetting=extended): ?

  • Yes, it works perfectly well in Maple 2016


Thanks to both of you, sorry for the inconvenience.


In my previous reply, please change


MAX := max(rhs~(Tally(SampleB)));
B0  := eval(B(rhs~(NE)[]), x=0);
scale := MAX/B0;





MODE := sort(Tally(SampleB), key=(x->rhs(x)))[-1];  #return the mode
XMAX := lhs(MODE);
HMAX := rhs(MODE):

B0  := eval(B(rhs~(NE)[]), x=XMAX);
scale := HMAX/B0;

This will give a better scaling of the graph of the probability function if the mode is not at x=0.
(try for instance r=2, b=1, d=0.1)


Last but not least: if the value of Xmax is large enough for the probability that
x > Xmax is very small, then the scale is simply given by

scale := `number of samples`

PS: you can verify that the sampling method described in the previous reply gives as good results as Statistics:-Sampling in the case d=0 (negative binomial distribution).



A spreviously quoted by Carl the case d=0 is obvious.
More precisely it corresponds to a negative binomial distribution with parameters r and b/(2+b)
(in Maple  NegativeBinomial(r, b/(2+b))  )

Carl also pointed the meaningless of the case d<0

I describe here the way to sample your  "pseudo distribution" B (which does not sum to 1 and then need to me normalized).
The last plot superimposes the probability function of B for r=3, b=2, d=1 (red line) and the histogram in discrete form.
I considered that only integers values of x have a sense.

For other valuers of r, b, d you will have to verify that the value of the variable Xmax is large enough (in this case the variable ShouldBe1 must have reached a "converged" value; to verify this increase the value of Xmax and look to what happened to ShouldBe1).

Let me know ih this suits you




A := (x+r-1)!/(r-1)!/x! * p^r * (1-p)^x



# verify A is the probability function of a NegativeBinomial distribution of parameters (r, p)

ProbabilityFunction(NegativeBinomial(r, p), x) assuming x > 0:
convert(%, factorial);



A := unapply(A, p):

pi := (b/2) / (1+d*x+b/2);
ip := (1+d*x) / (1+d*x+b/2);





# verify ip = 1-pi







B := unapply(A(pi) / (1+d*x), [r, b, d]);

proc (r, b, d) options operator, arrow; factorial(x+r-1)*((1/2)*b/(1+d*x+(1/2)*b))^r*(1-(1/2)*b/(1+d*x+(1/2)*b))^x/(factorial(r-1)*factorial(x)*(d*x+1)) end proc


Case d = 0

# When d=0 B(r, b, 0) is the probability function of a NegativeBinomial distribution of parameters (r, pi)
# where pi=b/(2+b)

simplify(subs(d=0, pi));

B(r, b, 0);
map(simplify, B(r, b, 0));







# sampling B(x, b, 0) is then trivial:
#  1/ set the values of r and b, for instance 3 and 5
#  2/ use Sample(NegativeBinomial(3, 5), 10)

Case d > 0
Note that the case d < 0 has no  meaning (see lso Carl's answers)

# numerical example (NE)

NE := [r=3, b=2, d=1]:



plot(B(rhs~(NE)[]), x=1..10)


# sequence of values of B(...) for integer values of x in [0, 100]

Xmax := 100:

PointwiseB := [ seq([x, B(rhs~(NE)[])], x in [$0..Xmax]) ]:

# given the plot above, assume the following quantity should be equal to 1
# if B was a true probability function

ShouldBe1 := add(op~(2, PointwiseB)):



# Normalize B:

PointwiseB := [ seq([x, B(rhs~(NE)[]) / ShouldBe1], x in [$0..Xmax]) ]:
evalf(add(op~(2, PointwiseB)));



# Draw the approximate cumulative function

c := CumulativeSum(op~(2, PointwiseB));

ApproximateCDF := zip((u, v) -> [u[1], v], PointwiseB, convert(c, list)):

plot(ApproximateCDF);    # not an exact representation... should be a stepwise function

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



# Sample ApproximateCDF:
#  1/ draw u from a Uniform(0, 1)
#  2/ find the

U := RandomVariable(Uniform(0, 1)):
ListTools:-BinaryPlace(op~(2, ApproximateCDF), Sample(U, 1)[1]);



SampleB := map(u -> ListTools:-BinaryPlace(op~(2, ApproximateCDF), u), Sample(U, 1000));

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



[0 = 550, 1 = 177, 2 = 92, 3 = 42, 4 = 24, 5 = 16, 6 = 13, 7 = 10, 9 = 5, 8 = 9, 11 = 7, 10 = 7, 13 = 3, 12 = 3, 15 = 3, 14 = 2, 19 = 1, 17 = 2, 22 = 4, 23 = 1, 20 = 1, 21 = 1, 27 = 1, 26 = 1, 31 = 3, 30 = 2, 36 = 1, 39 = 2, 33 = 3, 34 = 1, 44 = 2, 41 = 1, 40 = 2, 42 = 1, 63 = 1, 56 = 1, 73 = 1, 77 = 1, 67 = 1, 90 = 1, 94 = 1]


MAX := max(rhs~(Tally(SampleB)));
B0  := eval(B(rhs~(NE)[]), x=0);
scale := MAX/B0;

window := [0..20, default]:

      Histogram(SampleB, frequencyscale=absolute, discrete, thickness=10),
      plot(scale*B(rhs~(NE)[]), x=0..1, color=red, thickness=3),
      plot(scale*B(rhs~(NE)[]), x=1..100, color=red, thickness=3),








plot(scale*B(rhs~(NE)[]), x=0..1, adaptive=true)









Sorry, you're right, I forgot the "1/3" in the coefficient of x^3.

4 5 6 7 8 9 10 Last Page 6 of 17