mmcdara

4009 Reputation

17 Badges

6 years, 176 days

MaplePrimes Activity


These are answers submitted by mmcdara



A squared matrix which contains only one "1" on each line and each column and "0' everywhere else, is named a LHD matrix (LHD = Latin Hypercube Design).
I will say "LHD of 1" for reasons that will appear later.
Let N the dimension of such a matrix.

There are obviously N!  "LHD of 1", thus 6 "LHD of 1" in your case.

Now I'm going to replace, for one arbitrary "LHD of 1" among these six, some "0" by "1"

  1. Change nothing... 1 possibility (meaning "keep the original LHD")
     
  2. Build a single row with 2 "1" by starting f from the  "LHD of 1"
    This row can be chosen in 3 different ways and the "1" can be put in 2 different places.
    This lead to 3*2 = 6 different possibilities for a given "LHD of 1".
     
  3. Build two rows 2, both with 2 "1" by starting from the  "LHD of 1".
    These two rows can be chosen in 3 different ways (rows (1, 2), (1, 3) or (2, 3)).
    On the first row one can put the second "1" at 2 different places. Let c the number of the column this "1" is placed.
    There are two cases to consider:
    1. The second chosen row contains a "1" in column c: then there are 2 places left to place a second "1" on this row.
    2. The second row contains a "0" on column c: if one puts a "1" on column c, this column will necessarily contains 3 "1".
      So there is only 1 place left on the second row where a "1" can be placed.
    3. Finally there are 3*(2+1) = 9 possibilities to build two rows with 2 "1" per  "LHD of 1".
       
  4. Finally, build three rows with 2 "1" by starting from the  "LHD of 1".
    If one observe that this gicves a "LHD of 0", there are only 1 way to do that for a given "LHD of 1".
     

If we sum up these results, there are  6 * (1 + 6 + 9 + 1) = 102  matrices that meet your requirements


As they are probably different interpretations of your question, here is one
An_Idea.mw

 

The trick I used here is the transformation of each "Rieman" polygon into a prism by using plottools:-extrude.

Another method could be to construct  the piecewise function which represents the envelope of these "Rieman" polygons and make it turn (of course the rotating figure would then be a flat figure).



I think that exact integration is possible only if y = 0

Look at this

restart:
J := (y, tau) -> -2/3 + int(sqrt(x)/(exp((x-y)/tau ) + 1), x=0..infinity):
J(0, t) assuming t > 0;
J(0, t) assuming t < 0;
J(1, t) assuming t > 0;


In case you would be interested in numerical values only, you can do this

J := (y, tau) -> evalf(-2/3 + Int(sqrt(x)/(exp((x-y)/tau ) + 1), x=0..infinity, method=_d01amc)):
J(1, 1)
                          0.7297086143

 


1) use alias

restart:
alias(pode = plots:-odeplot):
p:= dsolve({D(y)(x) = y(x), y(0)=1}, type=numeric, range=-5..2):
pode(p);

2) use macro

restart
macro(pode = plots:-odeplot):
p:= dsolve({D(y)(x) = y(x), y(0)=1}, type=numeric, range=-5..2):
pode(p);


Wait for more skilled people to explain you the pros and cons of alias versus macro


I don't understand why "normality" as something to do with that?
Let A the data (a list or a vector).

This gives the interval of variation and the range

(min..max)(A);
(max-min)(A);

Instead of the last command you can also use this one

Statistics:-Range(A)

If you want more specific things like inter quantile range or quantile ranges, look to these terms in the Statistics package.
InterQuantileRange, Quantile, Range are not dependent on the underlying statistical distribution, so normality doesn't matter.

I seriously doubt that your Analytical_sol is the solution of [ode, bc].

Here is the correct formal solution:
(LargeExpressions is used ot ease Maple to calculate the solution)

with(LargeExpressions):
isolate(ode, diff(theta(y), y$2)):
rhs(%):
compact := collect(%, y, Veil[C]);
CS := [ seq( C[i] = Unveil[C]( C[i] ), i=1..LastUsed[C] ) ]:

                12                   11   3779136       10
 68024448 C[1] y   + 136048896 C[2] y   - ------- C[3] y  
                                             5            

      1259712       9   157464       8   209952       7
    - ------- C[4] y  + ------ C[5] y  + ------ C[6] y 
         5                25               25          

      5832       6   5832       5    81        4   27         3
    - ---- C[7] y  - ---- C[8] y  + ---- C[9] y  + --- C[10] y 
      125            125            1250           125         

       9         2   3            1       
    + --- C[11] y  + -- C[12] y + -- C[13]
      100            20           32      


dsolve([diff(theta(y), y$2)=compact, bc]):
Compact_Exact_Solution := rhs(%): 

# The solution is a 14 th degree polynomial in y degree(Compact_Exact_Solution, y) 
# while your Analytical_sol is a 4th degree polynomial in y

degree(Compact_Exact_Solution, y);

                               14


I would like to point out that if you had read the answer I had given to your initial question (which you destroyed without informing me), you would have had the (correct) answer to this one too.

To_Read_before_throwing_away.mw



linestyle is not correctly defined, proceed as you do for color :
 

colors := [black, red, blue, green, orange, cyan, purple, yellow, navy]; 
styles := [solid, longdash, spacedash, dash, dashdot]; 

for j to numelems(k__vals) do 
sol := dsolve(
         eval(
           [ode, bc], 
           [k = k__vals[j], We = .1, Br = .3, x = 0, Q = Q__vals[j], lambda = lambda__vals[j]]
         ), 
         numeric
       ); 
plt[j] := plots:-odeplot(
            sol
            , [y, theta(y)]
​​​​​​​            , y = 0 .. 1 
​​​​​​​            , axes = boxed 
            , color = colors[j] 
          # , style = [line] 
            , linestyle = styles[j] 
          # , symbol = [solidcircle, asterisk] # useless when style=line
          ) 
end do; 
plots:-display([seq(plt[j], j = 1 .. 5)])



By the way, the ode can be integrated formally.
(with a little help to Maple's 2015 version)
Help_graph_style_mmcdara.mw



Customize theattached file as you want (font, line thickness, characteristic vertical values and so on)

 

restart;

with(plots):

y1:=-1*v^2+1/8*v^4;

-v^2+(1/8)*v^4

(1)

hx, hy := 0.05, 0.05:

mima := [min, max]([solve(diff(y1, v))]);

pg := plot([[mima[2], eval(y1, v=mima[2])]], style=point, symbol=solidcircle, symbolsize=20, color=blue),
      textplot([mima[2], eval(y1, v=mima[2])-hy, "G"], align={'below'})
      :

[-2, 2]

(2)

X := 3:
p := plot(
  y1,
  v = -X .. X,
  #tickmarks=[0,0],
  color="Niagara Blue",
  labels= ["x", ""],
  labelfont = [TIMES, ITALIC, 16],
  size=[600,600],
  scaling=unconstrained,
  view=[default, eval(y1, v=mima[1])*1.1..eval(y1, v=X)*1.1]
):

 

Y := -1:
EF := sort(evalf([{solve(y1=Y)}[]]))[3..4], ["E", "F"];
pef  := seq(
          textplot([EF[1][n]+hx*(2*n-3), Y, EF[2][n]], align={'above', `if`(n=1, 'left', 'right')}),
          n=1..2
        ),
        seq(
          plot([[EF[1][n], Y]], style=point, symbol=solidcircle, symbolsize=20, color=black),
          n=1..2
        ),
        plot([[EF[1][1], Y], [EF[1][2], Y]], thickness=2, color=black),
        textplot([(EF[1][1]+EF[1][2])/2, Y, typeset(w = Y)], align={'above'})
        :

[1.082392201, 2.613125930], ["E", "F"]

(3)

Y    := 0:
CMD  := sort(evalf([{solve(y1=Y)}[]])), ["C", "M", "D"];
pcmd := seq(
          textplot([CMD[1][n]+hx*`if`(n=1, -1, +1), Y, CMD[2][n]], align={'above', `if`(n=1, 'left', 'right')}),
          n=1..3
        ),
        seq(
          plot([[CMD[1][n], Y]], style=point, symbol=solidcircle, symbolsize=20, color=red),
          n=1..3
        ),
        plot([[CMD[1][1], Y], [CMD[1][3], Y]], thickness=2, color=red),
        textplot([(EF[1][1]+EF[1][2])/2, Y, typeset(w = Y)], align={'above'})
        :

[-2.828427124, 0., 2.828427124], ["C", "M", "D"]

(4)

Y   := +1:
AB  := select(is, sort(evalf([{solve(y1=Y)}[]])), real), ["A", "B"];
pab :=  seq(
          textplot([AB[1][n]+hx*`if`(n=1, +1, -1), Y, AB[2][n]], align={'above', `if`(n=1, 'right', 'left')}),
          n=1..2
        ),
        seq(
          plot([[AB[1][n], Y]], style=point, symbol=solidcircle, symbolsize=20, color=green),
          n=1..2
        ),
        plot([[AB[1][1], Y], [AB[1][2], Y]], thickness=2, color=green),
        textplot([(EF[1][1]+EF[1][2])/2, Y, typeset(w = Y)], align={'above'})
        :
display(pg, pab, pef, pcmd, p, gridlines=true);

[-2.983115735, 2.983115735], ["A", "B"]

 

 

 


 

Download pot_mmcdara.mw

 

From  https://mathworld.wolfram.com/SimpleGraph.html :
"A simple graph, also called a strict graph (Tutte 1998, p. 2), is an unweighted, undirected graph containing no graph loops or multiple edges (Gibbons 1985, p. 2; West 2000, p. 2; Bronshtein and Semendyayev 2004, p. 346). A simple graph may be either connected or disconnected.".

Refering to Tutte the notion of simple graph only applies to undirected graph.

In case you are not satisfied with RandomDigraph (like @dharr  I use Maple 2015 too) you can replace this command with these lines

restart:
with(LinearAlgebra):
with(GraphTheory):
A := RandomMatrix(10, 10, generator=0..1, density=0.2);
A := A - DiagonalMatrix(Diagonal(A)):  # set A[i, i]=0 to avoid loops

G := Graph(A):
DrawGraph(G);

In case you want to create an random undirected graph to match Tutte's definition,  add these lines:

G_undirected := Graph(Vertices(G), convert~(Edges(G), set)):
DrawGraph(G_undirected);

Itis very likely that the image you display is a 2D projection of a 3D surface.

Here is an example (I arbitrary choosed Q as the second independent variable, be free to select abnother one if you want to).

plot3d(
  subs(We = .1, n = .5, lambda = 0.95e-2, x = 0, k = .1, u) 
  , y = 0 .. 1 
  , Q = 0 .. .6 
​​​​​​​  , axes = boxed 
​​​​​​​  , labels = ["y", "Q", "u"]
​​​​​​​  , style = surfacewireframe
  , shading=zhue
​​​​​​​  , orientation = [0, 90, 90]
​​​​​​​)


Putting aside the issue which concerns the execution time...

Maple's DataFrame is not that far from R's dataframe. Maple's DataFrame has less capabilities than R's dataframe, but the "logics of manipulation" are quite similar. Of course syntaxes differ, but if you are familiar with R, using dataframes in Maple should be relatively simple.
Just keep in mind that not all R's capabilities are avaliable in Maple's DataFrame.

A main issue in using xls/xlsx files is that ExcelTools:-Import and ExcelTools:-Export are operations which are rather slow.
My advice is to use binary file as soon as the number of data is large.

For instance

series(F, t = infinity) assuming theta > 0, mu > 1

It's not impossible that series eturn a result with other assumptions about alpha and beta (?)

Let E some event whose the true but unknown probability of occurrence is noted p.
A classical problem is the assessment of p.

Frequentist way
One realizes N trials of the phenomenon which can make E to appear.
Let K the number of times E has been observed.
The (basic) frequentist estimation of p is f=K/N.

Bayesian way
One still realizes N trials and observe K occurrences of E.
Obviously the value of K is a realization of some Binomial random variable Q with parameters N and q (q being any value between 0 and 1). It's clear that some values of q are more likely than others.
The bayesian way relies upon the knowledge of the values Q might take, and represents this knowledge in termes of a specific statistic distribution named Prior.
For Objective Bayesians, the classical choice is Prior(Q) = Beta(a, b) where a and b have, at least in the non hierarhical bayesian approach, fixed values.

One can show that a Beta distribution is a conjugate prior of a Binomial likelihood, which means that the posterior distribution of Q belongs to the same family than the prior.
Thus the main result : Post(Q) = Beta(a+K, b+N-k)

See for instance https://en.wikipedia.org/wiki/Conjugate_prior

So there is no need to use a Metropolis Algorithm here.
With  a=21, b=1, N=4, K=2Post(Q) = Beta(23, 3)

restart:
UseHardwareFloats := false
                             false
with(Statistics):
Mean('Beta'(23, 3));
evalf(%);
                               23
                               --
                               26
                          0.8846153846
Quantile('Beta'(23, 3), 0.025, numeric); 
                          0.7396941579
Quantile('Beta'(23, 3), 0.975, numeric);
                          0.9745346034

Remark: The convergence of your Markov Chain is slow because the your proposal (a Uniform(0, 1)) is highly inadequate.

I don't have time right now to go further but I bet that the posterior can be constructed in a thousand steps if the proposal is correctly chosen (as your prior concentrates in the range 0.8..1, take for instance a Uniform(-0.2/4, +0.2/4)). 

WATCHOUT : The prior being bounded, it is sure that some "prop" values you generate do not belong to [0, 1] (more certainly exceed 1). The expectation of the prior is 21/22 ~ 0.9545454545, meaning there are roughly 94% of chances that prop exceeds 1... look to this file
Download Beta+Uniform.mw

When you construct a Markov Chain an important quastion is "Did my chain behave correctly?"
A good looking chain must explore the space of possible values, must not be trapped on the same value for too many steps, its sample (after the burning phase) must posses a "correct" autocorrelation funcion and so on.
Just a question: have you observed the behaviour of your chain.

First: as you do not know a priori the number of iterations, it will be diffficult to predefine the array you talk about.
Secondly: the number of iterations depend on the method you use and on its parameters.
Thirdly: some methods have adaptive x steps and there is no simple relation between the number of the iteration and the physical time.

If yf1 and f2 do depend on the current iteration, why not write f1(x) and f2(x) instead ?

About events, look my answer to your previous question.

"Basically there is a very simply system of two differential equations and I would simply like to see an event triggered when the dependent variable is above 2"

In your worksheet you have written 

events = [[[0, 2-x < 0], [i(x) = i(x)+1]]]

but x is the independent variable.

Here are a few examples: the first and the second one correspond to what I understood (I arbitrarily used y(x) as the dependent variable used to monitor i(x)).

The third example is more complex: i(x) is augmented or reduce by one each time y(x) takes values +2 or -2. This example is purely illustrative.

The last example is a workaround to the condition "x > 2": it uses an auxuliary dependent variable w(x) defined by 

diff(w(x), x)=1, w(0)=0

Then, obviously, w(x)=x: "your" condition 

events = [[[0, 2-x < 0], [i(x) = i(x)+1]]]

then writes

events = [ [w(x)=2, i(x) = i(x)+1] ]

dsolve_events_mmcdara.mw

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