mmcdara

5635 Reputation

17 Badges

7 years, 189 days

MaplePrimes Activity


These are replies submitted by mmcdara

@Carl Love 

If I may say so, you are making an error of interpretation in saying "Using formal Statistics language, the p-value is less than 2%".

The p-value is used in inferential statistics and related to hypothesis testing.
Had @Christopher2222 provided empirical data collected among years, you could have tested the hypothesis "Was the name choosing really random?"
But it didn't, so the use of the name p-value is inadequate.

@delvin 

You don't have the same version of Maple so it's likely that the simplifications are not done the same way and that displaying orders differ.
As an exercice you can try to rewrite the expressions you got to get mine's.

Fot me the diffferences between us are not a subject worth wasting time on.

@delvin 

Can't you just be a little be patient and wait 75 seconds before panicking

Boy_in_a_hurry.mw

@Carl Love 

Yes of course, I wrote the relation I used without really paying attention.
Thanks for the correction.

I've proposed another approach meanwhile, see here
 

@Christian Wolinski 
Doing something like that (unfortunately the display time is too important) could be interesting.
The displays you will get correspond to those I gave in my answer vor diferent values of x.

Digits := 16:
p1 := plots[intersectplot](
    map(evalf, surface(eq1, y = 0 .. 20, z = -20 .. 20, x = 0 .. 1)),
    map(evalf, surface(eq2, y = 0 .. 20, z = -20 .. 20, x = 0 .. 1)),
    maxlev = 5, 
    grid = [100, 200, 20],  
    thickness = 3,  
    transparency = 0.3 
  ):

                               
Impacts := proc(X)
  plots:-display(
    p1, 
    plot3d(X, y = 0 .. 20, z = -20.. 20, color=black, transparency=0, style=surface)
    , view=[default, default, X..X+0.01]
    , scaling=constrained
    , orientation=[-90, 0, 0]
  );
end proc:

Explore(Impacts(X), parameters=[X = 0. .. 1.] );

In case you would be interested in improving this code?

@Carl Love @Christopher2222

Great Carl !

If I'm not mistaken a I believe to remember that a trick this kind is used to compute the Discrete Fourier Transform of a signal whose ength is not a power of 2.

Nevertheless one can simulate a discrete uniform random variable with support {1, .., 120} with only 2.9 throws ( a deterministic value, not an expectation of the number of throws)?
I'm not talking here of sampling a 6 players team out of a 120 roster for the procedure described below generates 10 numbers in one row.

The algotithm is:

  1. Step 1
    1. Roll dice:
      1. if its output is odd set B := 0, 3 (see in the attached file)
      2. if it is even set B := 1, 3.
    2. Roll the dice again, let T its output
      1. Set B6[1] := [T-1, B].
         
  2. Step 2 (T will always denotes the output of the throw)
    1. for i from 2 to 10:
      1. Roll the dice
        1. if T <= 3  set B := T-1  else  set B := T-4
      2. Roll the dice again and set B := T-1, B
      3. Roll the dice again and set B[i] := [T-1, B]
         
  3.  Each B[i] (i=1..10) is interpreted as a representation of an integer between 0 and 119 included.
    Convert each B[i]  in base 10 and add 1 to the result.
     

Remark
the upper bound 10 in the loop of step 2 comes from the fact there are 9 times more cards in the range 0..107 ( [0, 0, 0] ..[5, 5, 0] in base 6) than in the range 108..119 (( [0, 0, 3] ..[5, 1, 3] in base 6).
I'm not really proud of that because it forces the size of the card sample whose numbers belong to the range 1..108 to be exactly 9 times larger than the size of the card sample whose numbers belong to the range 109..120 (your algorithm let this emerge from randomness).

The result of this "enforced uniform distribution" strategy is that one needs exactly 29 throws to sample 10 cards.

 

restart:

with(Statistics):

F := proc({verbose::boolean:=false})
  local dice := rand(1..6):
  local B6, B10, throw_1, throw_2, throw_3, k, m:
  B6 := [[0$3]$10]:

  #-------------------- step 1
  k := 0:
  throw_1  := dice():
  k        := k+1:
  B6[1][3] := 3:
  if throw_1::odd then
    B6[1][2] := 0:
  else
    B6[1][2] := 1:
  end if:
  B6[1][1] := dice() - 1:
  k := k+1:

  #-------------------- step 2
  for m from 2 to 10 do
    throw_1 := dice():
    k := k+1:
    if throw_1 <= 3 then
      B6[m][3] := throw_1 - 1:
      B6[m][2] := dice() - 1:
      k := k+1:
      B6[m][1] := dice() - 1:
      k := k+1:
    else
      B6[m][3] := throw_1 - 4:
      B6[m][2] := dice() - 1:
      k := k+1:
      B6[m][1] := dice() - 1:
      k := k+1:
    end if
  end do:


  #-------------------- Base 10 conversion and card numbers
  B10 := map(b6 -> convert(b6, base, 6, 10), B6);
  B10 := map(b10 -> 1+add(b10 *~ [1, 10, 100][1..nops(b10)]), B10);

  `<,>`(B10), ` number of throws for 10 cards` = k
end proc:

randomize(170080957325093):
F();

Vector(10, {(1) = 112, (2) = 44, (3) = 76, (4) = 35, (5) = 63, (6) = 84, (7) = 9, (8) = 2, (9) = 51, (10) = 107}), ` number of throws for 10 cards` = 29

(1)

M := `<,>`(seq(F()[1], i=1..10^5));  # 10^6 samples
(min, max, Mean)(M);
Histogram(M, discrete=true, size=[1200, 400]);

sort(Tally(M), key=(x -> lhs(x))):

M := Vector(4, {(1) = ` 1000000 x 1 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

1, 120, Vector[row](1, {(1) = 60.4893560000000})

 

 

DU := RandomVariable(DiscreteUniform(1, 120)):
S  := Sample(DU, 10^6):

# Watch out: the next display will generate a 24 Mbyte file when saved
# Here is the image you will get


WantToPlot := false:

if WantToPlot then
DocumentTools:-Tabulate(
  [
    QuantilePlot(
      M, S
      , title="QQplot"
      , labels=["Dice", "Discrete-Uniform"]
      , labeldirections=[default, vertical]
    ),
    ProbabilityPlot(
      M, DU
      , title="PPlot"
      , labels=["Dice", "Discrete-Uniform"]
      , labeldirections=[default, vertical]
      , style = line
    )
  ]
  , width=50
);
end if;

 


Download Enhanced_result.mw


An interesting question is "Could it be possible to avoid trying another selection "if the card has already been selected" as Carl said.

Note that numerically sampling without replacement, let's say 10 cards, out of 120 is very easy:

combinat:-randperm([$1..120])[1..10]


For Christopher's son cards are not numeric but real and a very simple process mimics what the above command does.
Let's suppose the 120 cards are disposed on a circle in any order, figures visible.
A selection process, let's say yours, selects the card number N1. Write this number on a paper and turn upside down the card number N1.
Run your process a second time:

  • if the number N2 is different from N1, write this number on the same paper and turn upside down the card number N2,
  • if N2 = N1 chose the first card right card N1 which has its figure visible.

So there is no need to try another selection if the card has already been selected.

@delvin 

Did you read my file and did you understand it?
Given your last question I'm not sure of that.

BTW, did you check my code as I asked?

I t seems there is something wrong, either in the text you reproduce or in your 123.mw worksheet or in mine (please check it).
See here (set  verbose := true;   to see dthe details)

restart

eq := a^3*b*T^2*(diff(V(xi), `$`(xi, 2)))+a^3*b*T*(diff(V(xi), xi))+3*a^2*b*V(xi)^2-(3*a*c+2*b*omega)*V(xi)

a^3*b*T^2*(diff(diff(V(xi), xi), xi))+a^3*b*T*(diff(V(xi), xi))+3*a^2*b*V(xi)^2-(3*a*c+2*b*omega)*V(xi)

(1)


The "constant V" case:

# The "constant V" case:
# Let K a non null constant

Cst := factor(eval(eq, V(xi)=K)):

# K being non nul by definition one keeps onlu

Cst := Cst/K:

# Then Cst is null if

isolate(Cst, omega):

# The special case displayed in the excerpt you reproduce corresponds to K=a/3
# (why this choice?).
map(factor, eval(%, K=a/3))

omega = (1/2)*a*(a^2*b-3*c)/b

(2)


The "non constant V" case:

verbose := false:

Guess := xi -> (p[0]+p[1]*xi+p[2]*xi^(2)+p[3]*xi^(3))/(q[0]+q[1]*xi+q[2]*xi^(2)+q[3]*xi^(3))

proc (xi) options operator, arrow; (p[0]+p[1]*xi+p[2]*xi^2+p[3]*xi^3)/(q[0]+q[1]*xi+q[2]*xi^2+q[3]*xi^3) end proc

(3)

# Plug this guess into eq

f := eval(eq, V=(xi -> Guess(xi))):

# Under which conditions f=0 for any value of xi?

ind := convert(indets(f) minus {xi}, list);
sol := solve(identity(f, xi), ind):

if verbose then print~(sol): end if:

[T, a, b, c, omega, p[0], p[1], p[2], p[3], q[0], q[1], q[2], q[3]]

(4)

# You want to discard all solutions such that a*b*c = 0

sol_ok := map(s -> if eval(a*b*c, s) <> 0 then s end if, sol):
if verbose then print~(sol_ok): end if:

# You see there exist here solutions which give Guess(xi) = 0.
# As we supposed here that Guess(xi) is not a constant, we must
# discard them too:
sol_ok := map(s -> if has(eval(Guess(xi), s), xi) then s end if, sol_ok):
if verbose then print~(sol_ok): end if:

# What are the 3 different expressions of Guess?

Guesss_ok := map(s -> (eval(Guess(xi), s)), sol_ok):

if verbose then print~(Guess_ok): end if:

# Simplify now these expressions
simplify~(Guesss_ok)

[p[1]/q[1], p[2]/q[2], p[3]/q[3]]

(5)

# As they are all constants, there is a contradiction with the premisse that
# we are in the "non constant V" case

 

Download 123_mmcdara.mw

@Aung 

Oh sorry, I didn't pay attention (yhe 2D input mode is very luring and I never use it).
In Maple  exis not the exponential of xbut a name eraised to the power x(do evalf(e1)to convince you).
If mean "exponential of x" you must write

exp(x)

Concerning NonlinearFit , your syntax is not good, look to the help page.
Its first argument is the predictive model, here a expression which predicts the value of true_strain depending on the 7 parameters G[1], G[2], G[3], nu[12], tau[1], tau[2], tau[3].

Finally here are the solutions provided by NonlinearFit and Minimize.
They are practically identical but the extremely high value of the residual sum of squares means  

  • either that your model is dramatically wrong,
  • or thts is credible but that the two fitting methods are trapped in a local minimum.

So, ONE MORE TIME, provide either a phisycal range for each of the 7 parameters, or a reliable guess for each of them.
Why are you so reluctant to do this? Is this material top secret?

Here are the results obtained by the two methods:
prony_mmcdara_1.mw

@Aung 

for Optimization:-Minimize is simply

solution := Optimization:-Minimize(obj, constraints)

Unfortunately this lead to the error message

Error, (in Optimization:-NLPSolve) no improved point could be found

This is a quite common error, one way to fix it is to provide a "good" initial point.

A remark, doing this

indets(obj, name);
     {e, G[1], G[2], G[3], nu[12], tau[1], tau[2], tau[3]}

means obj doesn't depend on 6 parameters but on 8 parameters.
Adding constraints on e and nu[12] could help.
For instance, assuming no indeterminate is strictly negative, gives an answer (of course not physically acceptable):

solution := Optimization:-Minimize(obj, assume=nonnegative)
 [                    15                                  
 [6.122553512717879 10  , [e = HFloat(7.991226280161059), 

   G[1] = HFloat(0.0), G[2] = HFloat(0.0), G[3] = HFloat(0.0), 

   nu[12] = HFloat(0.0), tau[1] = HFloat(0.9999999472090075), 

   tau[2] = HFloat(0.9999999472090075), 

                                       ]
   tau[3] = HFloat(0.9999999472090075)]]

More of this a constraint of the form G[1] > 0 is very loose, G is likely a transverse Young modulus or something like that?
I supposed you know the material you did experiences upon and have some idea about the values of its mechanical  coefficients (I don't know, smething like G[1] >= 105,  G[1] <= 107  would help too.

Last but not least, maybe your objective function will have to be rewritten in a different way, maybe with other indeterminates too, in order to get solution.

BTW mmcdara is my personal account, sand15 is my professional one I use when at the office.

@Christian Wolinski 

Thanks for having seen my error.
I will see what can be done once corrected, likely not much fom a formal way.

Thanks again

@Christopher2222 

128 would have nee a very fair number of cards because  you could have gotten a card tossing 7 times a fair coin (each outcomes divides the card stack in two stacks of equal size). A simple dichotomy process.

But 120 is quite fair too.
Using the same coin, it happens that at some point (the 4th toss), you must have to split the stack of cards numbered 112 ro 128 into two stacks of size 8: stack 112..120 and stack 121..128.
If you face this situation do not toss the coin and chose the stack 112..120 and do three more tosses to select the card.
With 7 or 6 tosses your card is selected.
If you do not have a coin but a dice, you can decide that an odd face number corresponds to a head and an even to a tail.

Here is an extension of this strategiy which select the card with only "in mean" 4.373 throws (one can do a little bit better)

restart

with(Statistics):

SelectCard := proc()
  local dice, throw, pack, card:

  # split the stack of cards into 6 packs of size 20

  dice  := rand(1..6):
  throw := dice():
  pack  := [$(1+(throw-1)*20)..throw*20];

  # split each pack into 5 packs of size 4
  # convention: if throw=6  throx again the dice until you get a number from 1 to 5

  throw := dice():
  while throw=6 do throw := dice(): end do:
  pack  := pack[[$(1+(throw-1)*4)..throw*4]];
  # the dichotomy stage

  throw := dice():
  if throw::odd then
    pack := pack[1..2]:
  else
    pack := pack[3..4]:
  end if:

  throw := dice():
  if throw::odd then
    card := pack[1]:
  else
    card := pack[2]:
  end if;

  return card:
end proc:

# The expection of the number of throws is
# 4*5/6 + 5*1/6 + 6*1/36 + ...

4*5/6 + sum((4+i)/6^i, i=1..+infinity);
evalf(%);

328/75

 

4.373333333

(1)

Cards := Vector(10^5, i -> SelectCard())

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

(2)

Histogram(Cards, discrete=true, size=[1200, 400])

 

 

Download splitting_procedure.mw

As Christmas is coming I suggest you offer a set of Dungeons & Dragons dice to your son, he will be capable to select a card in just two throws :-)

@acer 

I realized this morning I've been a bit blunt and I apologize for this.

I agree that conversations with several persons often generate misunderstandings, not to say "conflicts".
I recently left a three-way thread because two different approaches were developing and this would only confuse the OP.

Sorry for my inappropriate attitude

@acer 

" ...but not actually requested in the original Question"

It wouldn't be the first time that a reply/answer develops around the initial question.
I even thought that was quite a common attitude among some people here? That's why I'm surprised by your remark.

@acer 

Your line 

G := [seq(1+irem(add(i*f(),i=1..n),n),j=1..50000)]:

doesn't fit what (I understand) @Christopher2222 wants to achieve. All the more weighting the output of a dice is not throwing a dice.

If it is to do exotic things, one can easily generate a discrete random variable over any interval with a single dice (but here again this is notv what I believe, @Christopher2222 wants ... either my solution below is in some sense optimal)

N := 10^5:
P := 20:
S := Sample(DiscreteUniform(1, 6), P*N):
[seq]([k, add(S[k+P*(n-1)], n=1..N) /(N*P) ], k=1..P):

# This last line generates a list of P lists of the form [k, pk], where k=1..P # and pk is the frequency of occurrence of k.
# It is easy to demonstrate that pk converges to 1/(P-1) as N --> +oo.
# Thus the distribution of this "one dice process" converges to the 
# distribution of a dicrete uniform random variable over {1, ..., P}

In fact a single run of this process can generate quasi discrete uniform distribution for any other value P' of P, provided that P' << N.
Examples are provided in the attached file where the lines above can be slightly changed (N and P unchanged) to generate a quasi discrete uniform distribution over {1,.., 40} and {1,.., 80}.
OneDiceProcess.mw

1 2 3 4 5 6 7 Last Page 2 of 119