490 Reputation

11 Badges

8 years, 72 days

MaplePrimes Activity

These are replies submitted by sand15

It seems to me that you wanted to answer the question EB1000 asked some hours ago?
Unfortunately, your answer turned out to be a question!


Thanks Tom,
 I do not understand either why "evalr" seems necessary. In any case, thank you for keeping me from going in circles.


Thank you acer.

As often there are several ways to tackle a problem and looking for the "optimal" solution seems to be an endless story.

By the way, I used myself the "subs({0="H", 1="T"},...)" method to return a solution in terms of "H" and "T.
Unfortunately writting
UseHardwareFloats := true:
V := floor~(SignalProcessing:-GenerateUniform(numThrows, 0, 2))
does not return integers making the {0="H", 1="T"} inoperative.
I tried the substitution {0.="H", 1.="T"}  but it didn't work neither.

I understand that floor, round, ceil do not return integers when the argument is a harware floating point number.
Is there a simple way to obtain, for instance, the integer part of harware floating point number?

@Carl Love 


Thanks Carl for your very documented answer.
I think I will  try to implement your suggestions, Not that the problem of calculating the probabilities of each sum deserves it, but rather out of personal satisfaction (in "practical" situation one would probably use Central Limit Theorems and the normal approximation when the number of dice exceeds a fex dozens).

By the way, I used in Q3 the package SignalProcessing and, from my student years I remember that convolution is generally be done, for effiency, by using FFT (2 directe + one inverse).
I do not know if SignalProcessing:-Convolution proceeds this way, but if it is the case, do you think its limitations you pointed out also impact FFT?



Have you thought about initializing the seed of the pseudo random generator?

This is a classical trap: executing twice the same worksheet generates exactly the same pseudo random numbers.To avoid this use:

Seed := randomize()

in a separate block justr after your top "restart" line

@Carl Love  Thanks for your remark.

It's rater difficult to compare the results produced by Q3 with the true ones. I did that for 3 and 4 dice and that was all good.
Obviously I was overconfident my assuming it would be same for higher values.

I suppose your observation definitely discredits Q3?


Maybe I'm wrong in my interpretation?
I understant A is a random variable with integer outcomes 1..6 and that each outcome has a different weight B

Would that be okay with you?
A := [$1..6];                   # A = [1, 2, 3, 4, 5, 6]
B := A *~ 10;                 # B = [10, 20, 30, 40, 50, 60]
P := B /~ add(B);           # P = [1/21, 2/21, ....6/21]
X := RandomVariable(ProbabilityTable(P));
N := 100;
Histogram(Sample(X, N), discrete);

# also:
plot(CDF(X, t), t=min(A)..max(A), discont=true)

Remark: if A was another set of outcomes, for instance A:=[1, 2, 3, 5, 8, 13], the plot of the histogram would be wrong (abscissas would be 1, 2, ..6). In this case you should correct the histogram plot by writting
Histogram(Sample(X, N), discrete, tickmarks=[[seq(i=A[i], i=1..6)], default])

Let me know if my answer is off the mark.


I agree: it's always better to try a simple approach.

Nevertheless this doesn't mean we can avoid some preliminary reflection.
When you write that the "penalty" function (I prefer personnaly the term of "objective" function) is  sum( (xmodel-xexp)^2+(ymodel-yexp)^2), you make implicit unwritten hypothesis.
This writting, while usual in least squares method, relies on underlyng assumptions:

  1. The disagreement between the model predictions (x*(t), y*(t)) and the experimental points (x(t), y(t)) come from random additive measurement errors.
  2. These errors are homoskedastic
  3. The measurement errors on X and Y are iid (independent and identically distributed)
  4. There are gaussian
  5. There is no measurement error on t
  6. There is no measurement error on alpha neither on v__0
  7. ...I'm not even sure this list is exhaustive

Under these hypotheses, sum( (xmodel-xexp)^2+(ymodel-yexp)^2) represents, up to ad additive constant, the log-likelihood of the experimental data given
the predictions from the model.

If only one of the previous assumptions fails, this writting is false.

In the case (3) fails, it's easy to modify your penalty function in  sum( ((xmodel-xexp)/sX)^2+((ymodel-yexp)/sY)^2) where sX and sY are the standard deviations
of the measurement errors on X and Y.

If the homoscedasticity assumption fails the penalty function should be  sum( wX(tn)*(xmodel(tn)-xexp(tn))^2+wY(tn)*(ymodel(tn)-yexp(tn))^2) where  wX(tn) and
 wY(tn) are weights depending on t (the previous case being a particular situation where   wX(tn) = (1/sX)^2 and wY(tn) = (1/sY)^2 for each tn.

if (4) fails the writting is simply false: think of errors with uniform distributions!

if (5) fails the only correct way to account for measurement errors on t is a bayesian approach.
The often invoked "generalized least squares method" is not rigourous an would fail if (4) fails too.

I stop here.
Let me be clear, I'm not saying your work is meaningless. I just say that writtig the formula you used to express a "penalty" function has no sense if you
do not explicitely state the hypothesis that makes it reliable.
Too many people use statistical methods as if they were cooking recipes.
I thought this clarification was necessary in order Erik may form himself a correct idea about the true complexity of its problem, and thus about the most suited method to use in order to solve it

@Carl Love 

The OP wrote "pdf(P(x)=1/k+1  x=0..k) "

Does this mean "the pdf of a random variable (RV) X such that Prob(X=x) = 1/k + 1 for each value x in 0..k" ?
If it is so, then X is not a RV because, obviously the sum of all Prob(X=x) over 0..k doesn't sum to 1.

Does this mean "the pdf of a RV Xk is defined by the function P(x) = 1/k + 1 for each value x in 0..k" ?
If it is so, then X is not a RV because, here again the sum of all P(x) over 0..k doesn't sum to 1.

That is why I interpreted "pdf(P(x)=1/k+1  x=0..k) " the way I did in my reply.
I do not pretend I was right, but the original question is not correctly formulated and leaves open many doors.

Following a suggestion from acer about one of my previous posts concerning the PDF of the sum of two RVs (usage of 'inert') I tried this

roll := rand(0. .. 1.):
z := (a, b, c) -> Mean(c*RandomVariable(Uniform(a, b))):
y := (a, b, c) -> evalf(Mean(c*RandomVariable(Uniform(a, b)), inert)):

abc := roll(), roll()+1., roll:

z(abc);  # fail almost surely with probability 1
y(abc);  # always returns the good result

I remember acer writing he was going to submit a bug report, maybe the fact z(abc) fails could be related to the same bug?

@Carl Love 

This last file contains:

  1. A procedure named 'MD' which does basically the same thing that the procedure 'md' did in my previous sending, but which uses a more formal approach.
    The plot presents no cusp, which probably means there remains an error in my 'md' procedure
  2. An second approach based on a polar representation of the ellipse where the origin of the radius vector rho starts from the origin and not from one of the focal points.
    I use here the following definition : MeanDistance := int(rho(theta), theta=alpha, beta) / (beta-alpha)

The last plot compares the two definitions of the mean distance to the origin.

If I had to choose one of them I would prefer the second one.

Not is all rosy, for this second computation leads to some numerical problems; in particular the formal solution is undefined for alpha (beta, theta) = 0, Pi, 2*Pi

@Carl Love 

No, I can't

I was myself dubious about this cusp: after all one computes an integral and a cusp would mean that the integrand is not continuous (which it is not).

I wrote in my previous answer "(if I'm not mistaken)", so I came back to my code 'md' and found a mistake for the case "alpha > Pi" (==> beta > Pi).
I found this by comparing the plots for

  1. alpha = Pi/2 -epsilon, beta=Pi/2+epsilon
  2. alpha  = 3*Pi/2 -epsilon, beta=3*Pi/2+epsilon

After correction these two plots are the same up to a translation of Pi.

Nevertheless the cusp still remains when I plot the mean distance between alpha = 0 and beta in 0..2*Pi

This looks all the most strange that the same plot for alpha = 3*Pi/2 -epsilon and beta in 3*Pi/2-epsilon,..3*Pi/2+epsilon 
doesn't reveal any cusp.
I suspected it could have been a problem related to the value of 'numpoints', but it is definitively a wrong track

I propose you to investigate this point deeper when I'll be home (in about 8 hours).
In the meantime here is the new procedure.

@Preben Alsholm 

Thanks Preben



@Carl Love

Thank you all for all for all your contributions.

I also undertood that the format of the files is not neutral, nor even their order of creation

The first reply from vv made me questioned myself about the results returned by addressof.
I didn't know of this command and I did a few things around my initial issue.
They are summarized in the attached file


  1. To refer to the last vv's reply...
    Yes the addresses of A, B, and C are different (please look to the code below).
    But as B:=eval(A) makes a strong link between A and B, I could have thought the their adresses would have be the same (contrary to those of A and C:=copy(A)).
    The fact that changing one ientry of A (B) modifies accordingly B (A) seemed to advocate this intuition.
    I was definitively wrong

    PS: for those of you who know FORTRAN; there was a statement named EQUIVALENCE, which enabled to access the same memory location through different names
  2. Strangely for me, the adresses of A[n], B[n] and C[n] are the same for each index n (please look the code below)
    Why this?
    Could it be that addressof(A[1]) does not return the address of A[1] as I was expected for?
  3. Plus a few auxiliary quesions  you will discover by your own

If you get a little bit of time, could you be kind enough for enlighten me about the few points I raise here
Nothing critical for me, just the desire to understand a little bit more deeply how MAPLE works.




A := table([1="a", 2="b"]):
B := eval(A):
C := copy(A):

# All addresses are different
# A and C being different objects, I understand why the address of C is not the one of A.
# But I thought that B:=eval(A) was like some kind of "alternate naming of A" and that the
# addresses of A and B were the same (which seemed consistent with my changing B accordingly
# changed A).
# I was obviously wrong.

addressof(A) , pointto(addressof(eval(A)));
addressof(B) , pointto(addressof(eval(B)));
addressof(C) , pointto(addressof(eval(C)));


18446744074328626430, table( [( 1 ) = "a", ( 2 ) = "b" ] )


18446744074328627326, table( [( 1 ) = "a", ( 2 ) = "b" ] )


18446744074328627422, table( [( 1 ) = "a", ( 2 ) = "b" ] )



# An auxiliary question...

addA := [disassemble(addressof(A))];   # what do these addresses mean ?
addB := [disassemble(addressof(B))];   # ""
addC := [disassemble(addressof(C))];   # ""


seq(print(addA[k], pointto(addA[k])), k=2..4);  # returns results I don't understand, what do they mean ?

[8, 18446744074328626590, 18446744074328154110, 18446744073709551745]


[8, 18446744074328626590, 18446744074328154110, 18446744073709551747]


[8, 18446744074328627902, 18446744074328154110, 18446744073709551749]



18446744074328626590, table( [( 1 ) = "a", ( 2 ) = "b" ] )




18446744073709551745, 65


# Adresses of the first index are the same despite the fact that addressof(A) = .. = addressof(C)
# Would have thought addressof(A[1]) = addressof(B[1]) <> addressof(C[1])

locA1ini := addressof(A[1]):  locA1ini, pointto(addressof(A[1]));
locB1ini := addressof(B[1]):  locB1ini, pointto(addressof(B[1]));
locC1ini := addressof(C[1]):  locC1ini, pointto(addressof(C[1]));

18446744074328626462, "a"


18446744074328626462, "a"


18446744074328626462, "a"


# Adresses of the second index are the same despite the fact that addressof(A) = .. = addressof(C)
# Would have thought addressof(A[2]) = addressof(B[2]) <> addressof(C[2])

addressof(A[2]) , pointto(addressof(A[2]));
addressof(B[2]) , pointto(addressof(B[2]));
addressof(C[2]) , pointto(addressof(C[2]));

18446744074328626494, "b"


18446744074328626494, "b"


18446744074328626494, "b"


# Make a change of B[1] and verify that A[1] is changed accordingly and C[1] remains unchanged

B[1] := "B[1] has changed":

table( [( 1 ) = "B[1] has changed", ( 2 ) = "b" ] )


table( [( 1 ) = "B[1] has changed", ( 2 ) = "b" ] )


table( [( 1 ) = "a", ( 2 ) = "b" ] )


# adresses of A[1] and B[1] have changed but they are still equal
# address of C[1] is kept unchanged
# A table being a mutable structure I would have thought the address of B[1] would have kept unchanged.
# Could it be due to the fact that the numbers of bytes of "a" and B[1] being different, a reallocation
# of the internal memory is necessary?

printf("%a %a %a\n", locA1ini, addressof(A[1]) , pointto(addressof(A[1])));
printf("%a %a %a\n", locB1ini, addressof(B[1]) , pointto(addressof(B[1])));
printf("%a %a %a\n", locC1ini, addressof(C[1]) , pointto(addressof(C[1])));

numelems(convert("a", bytes));
numelems(convert(B[1], bytes));

18446744074328626462 18446744074329309118 "B[1] has changed"
18446744074328626462 18446744074329309118 "B[1] has changed"
18446744074328626462 18446744074328626462 "a"







# Make a change of C[1] and verify that A[1] and B[1] are kept unchanged

C[1] := "C[1] has changed":

table( [( 1 ) = "B[1] has changed", ( 2 ) = "b" ] )


table( [( 1 ) = "B[1] has changed", ( 2 ) = "b" ] )


table( [( 1 ) = "C[1] has changed", ( 2 ) = "b" ] )


# adresses of A[1] and B[1] are kept unchanged
# only the address of C[1] has changed

printf("%a %a %a\n", locA1ini, addressof(A[1]) , pointto(addressof(A[1])));
printf("%a %a %a\n", locB1ini, addressof(B[1]) , pointto(addressof(B[1])));
printf("%a %a %a\n", locC1ini, addressof(C[1]) , pointto(addressof(C[1])));

18446744074328626462 18446744074329309118 "B[1] has changed"
18446744074328626462 18446744074329309118 "B[1] has changed"

18446744074328626462 18446744074329311198 "C[1] has changed"






Thanks Tom for your help.
Filtering the original matrix was not my primary intention but it could be useful in the future.



@Carl Love 

indets(a, specfunc(VIEW))[]; is indeed a more general way to get the information
Thank you


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