Advanced Search

mmcdara

7891 Reputation

22 Badges

9 years, 55 days
Contact mmcdara

MaplePrimes Activity


These are replies submitted by mmcdara

@Kitonum Think I should change my s...

mmcdara 7891
March 15 2019
0 0

@Kitonum 

Think I should change my spectacles.

Sorry for the disagreement

@acer @RimAgree with Acer.More of t...

mmcdara 7891
March 15 2019
0 0

@acer @Rim

Agree with Acer.

More of this the fonction f acting locally on a each element of M, the problem is perfectly scalable. So if you have N material nodes you can expect at most  to divide the computation time by N.
So if you expect very large speedups on a "classical"l 4 nodes machine, using Grid probably won't be a solution

By experience, if you work with Windows 8 or higher, the OS seems to have its own strategy to "charge" each node the way it wants. This seems to have a slight interference with Grid (I observed better speedups with Windows 7 than with Windows 8, but you can bet on a speedup of about N).
If you have hyperthreaded nodes, the speed up is between 1.3N and 1.6N.
 

Hi Tom, this is a problem I rather often...

mmcdara 7891
March 15 2019
0 0

Hi Tom, this is a problem I rather often come across.
For info, I work here with Maple 2015.2, under Mac OS Mojave 

Other info, I talked about this problem a few days ago (not about is specifically, just a PS at the end of some reply)
@acer Acer, I realized that this coding&nbs...

@tomleslie Your result (2) is alway...

mmcdara 7891
March 14 2019
0 0

@tomleslie 

You write "thus x__i is a single 'name' which is completely independent of the summation index 'i' . You might as well write this as the name 'xi'.", and I agree.

But as a consequence your result (2) is always equal to 0 whatever the values of N and xi

Don't you think the procedure g(N, k) below is closer to what Hammur has in mind ?


 

> 

restart:
> 

g := (N, k) -> N+k*(Sum(ln(x[i]-B), i = 1..N))-N*k*(Sum((x[i]-B)^k*ln(x[i]-B), i = 1..N))/(Sum((x[i]-B)^k, i = 1..N));

proc (N, k) options operator, arrow; N+k*(Sum(ln(x[i]-B), i = 1 .. N))-N*k*(Sum((x[i]-B)^k*ln(x[i]-B), i = 1 .. N))/(Sum((x[i]-B)^k, i = 1 .. N)) end proc

 (1)
> 

diff(g(N, k), k);

Sum(ln(x[i]-B), i = 1 .. N)-N*(Sum((x[i]-B)^k*ln(x[i]-B), i = 1 .. N))/(Sum((x[i]-B)^k, i = 1 .. N))-N*k*(Sum((x[i]-B)^k*ln(x[i]-B)^2, i = 1 .. N))/(Sum((x[i]-B)^k, i = 1 .. N))+N*k*(Sum((x[i]-B)^k*ln(x[i]-B), i = 1 .. N))^2/(Sum((x[i]-B)^k, i = 1 .. N))^2

 (2)
> 

simplify(evalf(diff(g(1, k), k)));

0.

 (3)
> 

simplify(evalf(diff(g(2, k), k)));

(-1.*ln(1.*x[1]-1.*B)*((x[1]-B)^k)^2+ln(1.*x[1]-1.*B)*((x[2]-B)^k)^2+ln(1.*x[2]-1.*B)*((x[1]-B)^k)^2-1.*ln(1.*x[2]-1.*B)*((x[2]-B)^k)^2-2.*k*ln(1.*x[2]-1.*B)^2*(x[2]-B)^k*(x[1]-B)^k-2.*k*ln(1.*x[1]-1.*B)^2*(x[1]-B)^k*(x[2]-B)^k+4.*k*(x[1]-B)^k*ln(1.*x[1]-1.*B)*(x[2]-B)^k*ln(1.*x[2]-1.*B))/((x[1]-B)^k+(x[2]-B)^k)^2

 (4)
> 

 


 

Download Correction.mw

@acer Acer, I realized that this co...

mmcdara 7891
March 12 2019
0 0

@acer 

Acer, I realized that this coding 

  samplingMethod := proc(m,s,N)
    local U, G:
    uses ST=Statistics;
    U := ST:-RandomVariable(Uniform(0, 1));
    G := ST:-RandomVariable(Normal(1/2, 1));
    m +~ s /~ ( ST:-Quantile~(G, ST:-Sample(U, N)) )^~2
  end proc;

is "better" than this one

  samplingMethod := proc(m,s,N)
    local U, G:
    uses ST=Statistics;
    m +~ s /~ ( ST:-Quantile~(':-Normal'(0, 1),
                             ST:-Sample(':-Uniform'(1/2, 1), N)) )^~2
  end proc;

When you draw a sample of size N this second coding generates N+2 random variables while the first one only generates 2.
More of this the first one is twice faster...

Download Improvement.mw



PS: I sometimes get this message when I try to upload the content of a worksheet. Might you have any idea why this happens?
Maple Worksheet - Error
Failed to load the worksheet /maplenet/convert/Improvement.mw .

 

 

 

@acer Sorry I minterpreted.Yes, the...

mmcdara 7891
March 11 2019
0 0

@acer 

Sorry I misinterpreted.
Yes, the support is open(m)..open(+infinity) and you did right.

The Mean (and the variance) are infinite (I just added ':-Mean' = +infinity.

Effectively Quantile~(...) is a slow operation. I often use this "trick" *** to sample a truncated Normal distribution, option 'numeric' helps reducing the computation time but the operation remains costly.
Using envelope in the sampling operation doesn't seem to be a good workaround for sampling truncated distribution (as a rule the support of the envelope should contain the one of the target distribution).

By the way, I didn't go deep to find out how to sample a Levy's law: I just translated word for word the code used in R. For the moment I'm working on Stable distributions whichthe Levy one is a particular case. I will compare the sampling time of a Levy distribution with the one of the corresponding Stable distribution (which makes no use of Quantile).

Thanks again


*** generate a sample U from Uniform(p[1], p[2], then apply Quantile~(Normal(0, 1), U) to generate a sample in a range [q[1],q[n]]  where q[n] is such that Probability(G < q[n]) = p[bn] ig G ~ Normal(0, 1)

@acer Definitely excellent !Support...

mmcdara 7891
March 11 2019
0 0

@acer 


Definitely excellent !

Support is wrong but I added  ':-Support'= m..+infinity' to fix the problem.
No "funny business" happens with option remember kept in place.

Thanks acer


PS: radaar asked a few days ago for a Levy distribution sampler.
I gave him an answer but, once again, I faced this same inability  to build correctly a RV with non instanciated parameters. I then took this occasion to open this thread here.
I'm forwarding your work to radaar, it might interest him

@Teep  Given tomleslie's and v...

mmcdara 7891
March 10 2019
0 0

@Teep 

Given tomleslie's and vv's remark, here is the adptation of Graph_2.mw to a DIRECTED graph (Graph_2.mw adresses onlye the cas of an UNDIRECTED graph)

Graph_3.mw

@vv  tomleslie made me remark that...

mmcdara 7891
March 10 2019
0 0

@vv 

tomleslie made me remark that the matrix A wasn't symmetric, so you were" right, the graph is probably directed.
I was lured by the term "connection" which (to my opnion) seemd to refer to a "symetric connection").
I would have been more careful...

@Teep  Assuming you have an UNDIRE...

mmcdara 7891
March 10 2019
0 0

@Teep 

Assuming you have an UNDIRECT graph (contrary to vv):

As I mention in the worksheet the rendering depends on the version of Maple you use (for instance Maple 2018 returns the pictures vv displayed)

> 

restart
> 

with(GraphTheory):
> 

A := Matrix(4, 10, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 1, (2, 5) = 0, (2, 6) = 1, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 1, (3, 1) = 1, (3, 2) = 1, (3, 3) = 0, (3, 4) = 0, (3, 5) = 1, (3, 6) = 0, (3, 7) = 1, (3, 8) = 1, (3, 9) = 1, (3, 10) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 0}):
> 

edges := map(x -> if A[op(x)] = 1 then {x[]} end if, {indices(A)}):
G := Graph(edges):
DrawGraph(G, style=spring);

 

You will have a different rendering depending on the version you use, here is mine

> 

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

 (1)
> 

LocalConnections := proc(G::Graph, n::integer, col1::symbol:=red, col2::symbol:=gray, how::symbol:=spring)

  # default colors are red and gray, default style=spring

  local g, ie:

  if member(n, Vertices(G)) then
    g  := CopyGraph(G):
    ie := IncidentEdges(G, n):
    HighlightEdges(g, ie, col1):
    HighlightEdges(g, Edges(g) minus ie, col2):
    print(DrawGraph(g, style=how));
  else
    error cat(n, " is not a vertex of the input graph");
  end if:

end proc:
> 

# a few different uses
#
# LocalConnections(G, 2, red, gray, spring)
# LocalConnections(G, 2, red, gray)
  LocalConnections(G, 2);


# You can keep customizing this procedure by coloring vertices:
#
#   HighlightVertex(g, n, SomeColor):
#   HighlightVertex(g, Neighbors(g, n), AnotherColor):
#   HighlightVertex(g, [({Vertices(g)[]} minus {Neighbors(g, n)[]} minus {n})[]], StillAnotherColor): 

 

 
> 

PartialLocalConnections := proc(G::Graph, n::integer, L::{list, set}, col1::symbol:=red, col2::symbol:=gray, how::symbol:=spring)

  # default colors are red and gray, default style=spring

  local g, ie, nn:

  if member(n, Vertices(G)) then
    g  := CopyGraph(G):
    ie := {}:
    for nn in L do
      if member({n, nn}, Edges(g)) then
        ie := ie union {{n, nn}};
      else
        WARNING(cat("There is no ", {n, nn}, "edge in the input graph"));
      end if:
    end do:
    if ie <> {} then
      HighlightEdges(g, ie, col1):
      HighlightEdges(g, Edges(g) minus ie, col2):
      print(DrawGraph(g, style=how));
    end if;
  else
    error cat(n, " is not a vertex of the input graph");
  end if:

end proc:
> 

PartialLocalConnections(G, 2, [3, 4, 6, 9])

Warning, There is no {2, 9}edge in the input graph

  

 
> 

 


 

Download Graph_2.mw

 

@vv  Why do you think the graph sh...

mmcdara 7891
March 10 2019
0 0

@vv 

Why do you think the graph should be directed?

@Carl Love  Thanks Carl, works per...

mmcdara 7891
March 09 2019
0 0

@Carl Love 

Thanks Carl, works perfectly.

I thumb up

@Carl Love  Nice. Carl, I'm p...

mmcdara 7891
March 09 2019
0 0

@Carl Love 

Nice.

Carl, I'm presently using Maple 2015 (so no Iterator package and then no hope to run your code... I'll do that later with Maple 2018), but I'm confused with the line 

P[(fX:= f(seq(X)))]:= P[fX] + mul(Pr(R[r]=X[r]), r= 1..nR)

Furthermore the first ':='generates an < ':=' unecpexted > error
Could you explain me what this line does and how it would write in Maple 2015?

Tnaks in advance

@erik10 @acerTo complete my previou...

mmcdara 7891
March 07 2019
0 0

@erik10 @acer


To complete my previous post just type this

M := outcome -> Matrix(6, 6, (i,j) -> `if`(g(i,j)=outcome, 1, 0)):
MAT := Matrix(2, 3, [ seq(plots:-matrixplot(M(n), heights=histogram, title=n), n in outcomes) ]):
DocumentTools:-Tabulate(MAT, width=60)

 


What I integrate are the "2d histogram" functions matrixplot exhibits

 

A workaround to compute probabilities...

mmcdara 7891
March 07 2019
0 0

@acer  @erik10

Let f(A, B) the function of the 2 RVs A and B (I used "my" Heaaviside form but I guess acers' K function could be handled too)

One can can compute the probabilities by explicitely doing a double integral of the characteristic function of f(A, B)=n for al possible outcomes of f (-4, 0, 1, 2, 3, 4)

I use here Student[MultivariateCalculus][ApproximateInt] in such a way that the integral is exact (partition=[6, 6], method=midpoint, plus an adhoc definition of the upper and lower integration bounds.

I did not encapsulate all the stuff within some "Probability" procedure in order to keep the things readable.


 

> 

restart:
> 

with(Statistics):
> 

Dice1 := RandomVariable(DiscreteUniform(1, 6));
Dice2 := RandomVariable(DiscreteUniform(1, 6));

_R

  

_R0

 (1)
> 

g := (Dice1, Dice2) -> -4*(1-Heaviside(Dice1-3/2)*Heaviside(Dice2-3/2)) + abs(Dice2-Dice1)*Heaviside(Dice1-3/2)*Heaviside(Dice2-3/2)

proc (Dice1, Dice2) options operator, arrow; -4+4*Heaviside(Dice1-3/2)*Heaviside(Dice2-3/2)+abs(Dice2-Dice1)*Heaviside(Dice1-3/2)*Heaviside(Dice2-3/2) end proc

 (2)
> 

with(Student[MultivariateCalculus]):
> 

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

[-4, 0, 1, 2, 3, 4]

 (3)
> 

UpperBound := outcomes +~ 0.5;
LowerBound := UpperBound -~ 1;

[-3.5, .5, 1.5, 2.5, 3.5, 4.5]

  

[-4.5, -.5, .5, 1.5, 2.5, 3.5]

 (4)
> 

pdf := NULL:

for n from 1 to 6 do
   chi := Heaviside( g(u,v) - LowerBound[n] ) - Heaviside( g(u,v) - UpperBound[n] ):
   ApproximateInt(chi, u=0.5 .. 6.5, v=0.5 .. 6.5, method=midpoint, partition=[6, 6]):
   pdf := pdf, %;
end do:

proba := convert~([pdf]/add([pdf]), rational);

[11/36, 5/36, 2/9, 1/6, 1/9, 1/18]

 (5)
> 

ExactProba := [11/36, 5/36, 2/9, 1/6, 1/9, 1/18]

[11/36, 5/36, 2/9, 1/6, 1/9, 1/18]

 (6)
> 

 


 

Download Probabilities.mw


 

First 131 132 133 134 135 136 137 Last Page 133 of 154