sand15

520 Reputation

10 Badges

4 years, 355 days

MaplePrimes Activity


These are replies submitted by sand15

@radaar 

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?
restart:
with(Statistics):
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.

@tomleslie 

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  https://en.wikipedia.org/wiki/Homoscedasticity
  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

restart:
with(Statistics):
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

Ellipse3b.mw

@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.

Ellipse2.mw

@Preben Alsholm 

Thanks Preben

@vv 

@tomleslie 

@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

Basically:

  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.

 


 

restart:

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

print();
 

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

 

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

 

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

 

(1)

# An auxiliary question...

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

print();

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" ] )

 

18446744074328154110

 

18446744073709551745, 65

(2)

# 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"

(3)

# 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"

(4)

# 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":
eval(A);
eval(B);
eval(C);

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

 

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

 

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

(5)

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

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

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

 

 

1

 

16

(6)

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

C[1] := "C[1] has changed":
eval(A);
eval(B);
eval(C);

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

 

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

 

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

(7)

# 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"

 

 


 

Download AboutTables.mw

@tomleslie 

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

@Kitonum 
Thanks

 

@Carl Love 

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

 

@Carl Love 

This is a workaround. Do you know what that means? It is not a denial of the existence of a bug! It acknowledges the bug's existence and provides the user with an easy and practical way to "work around" the bug until it is fixed.

Yes, I know what a workaround is. It can be used in any of these two situations:

  1. You face a Maple impossibility to do a specific operation and you replace it by something else.
    Nothing can answer to all issues, every product, here a Maple, has its own limitations and using it beyond its frontiers may require workarounds.
     
  2. You face a bad behaviour of Maple, let's say a "bug", even if it is a taboo word here, and you proceed otherwise to avoid the bug.
     

The main problem with the most part of workarounds published here (a significant part of all the replies), is that the op doesn't know if they are proposed to help using Maple out of its application domain, or to circumvent a bug (as a rule workarounds are not explained ... even if you are among those that constitute an exception).
Given the taboo character of this word, there seems to be a collective denial to their existence, which, for me,  is definitely unbearable.
For many years I develop computational codes and I face, as anybody does too, user feedbacks about bugs. I'm not ashamed at that.

This is why I'm happy to read from you  It is not a denial of the existence of a bug!
Not so often acknowledge here.

After my reply to the first vv's answer you wrote  I suspect that you didn't type his commands into your Maple 2018
What vv proposed was this
         restart:
         with(Statistics):
         N :=3; 
         X := RandomVariable(Binomial(N, 1/2)):
         F:=CDF(X, s):


It solved nothing at all as it is confirmed by the the attached file (or I did some mistake in this file I'm blind to find), so I'm not sure you were right in saying I suspect that you didn't type his commands into your Maple 2018

Major_Hassle.mw


Far from any controversy now.
But I think that Vv is right about you: You're more interested in arguing than getting practical information on how to use Maple
About my propensity of arguing ...  you can be right. But, on the other side, if the substance of the discussions could let me think that there was no denial to recognize the existence of bugs, I would be  probably less vindictive.
 

PS: Your memory is better than mine for I do not remember what my first post was.

@acer 

First of all: thank you for these "dispassionnate" answers.

I have a few questions ad remarks:

About your first reply:

  1. F := value(CDF(X, s, inert)) returns a more credible expression than CDF(X, s): do you intend to modify the CDF procedure in order that it will do this automatically?
    All the more that, if I'm not mistaken, there is no advice in the help pages to proceed the way you did.
  2. F := value(CDF(X, s, inert)) F := sum(piecewise(_t < 0, 0, binomial(3, _t)*(1/2)^_t*(1/2)^(3-_t)), _t = 0 .. s) is close to the Maple 2015 output of  CDF(X, s) piecewise(s < 0, 0, 3 <= s, 1, sum(binomial(3, i)*(1/2)^i*(1/2)^(3-i), i = 0 .. floor(s))).
    Nevertheless  some discrepancies remain :
    1. value(CDF(X, s, inert)) explicitely tells nothing about what happens beyond s>3 (contrary to Maple 2015).
    2. The upper summation bound is no longer floor(s) but s, even if it seems to be exactly the same thing, why this change?



About your second reply:

  1. F := CDF(X, s) assuming s <= 3 doesn't suit me for it implicitely tells that the CDF is not defined for values of s strictly greater than 3.
     
  2. Given the fact that PDF(X, t) is a linear combination of "Dirac functions", the computation of the cdf through the command   cdf := int(P(X, t), t=-infinity..s) provides a linear combination of Heaviside functions (exactly what Maple 2015 does too).
    So, why does this Psi function come from ?


 

Finally

  1. Do you think of recovering the old Maple 2015 CDF procedure ?
     
  2. For my sake, would you qualify as bug the fact  CDF(X, s) is not returning a proper expression?

 


CDF of Binomial(6, 1/2) is OK  but not CDF of Binomial(6, 1/4) 
The step at location [0, 1[  (or ]0, 1] because french and english customs do not agree on this point) is not present!

X := RandomVariable(Binomial(6, 1/4)):   # or Binomial(6, 0.25) if you prefer
F := CDF(X, s):
plot(F, s=-1..7, discont=true, gridlines=true, axis[1]=[gridlines=7]);


The CDF has an uncomprehensible form invoking the Psi function. 
Could anyone from the development team restore what was done in Maple 2015?

@acer 
I'm so sorry acer for not having replied sooner.
After this inexcusable delay I greatly thank you for all these precious informations.

2 3 4 5 6 7 8 Last Page 4 of 14