mmcdara

787 Reputation

12 Badges

3 years, 136 days

MaplePrimes Activity


These are answers submitted by mmcdara

Here are two ways to do this. 
 


 

restart:

with(Statistics):

C := RandomVariable(Cauchy(0, 1)):

# The easiest way

N := 10^4:
C_sample := Sample(C, N, method=[envelope, range=-1..1]):
Histogram(C_sample);

[Mean, Variance](C_sample);

CodeTools:-Usage( Sample(C, N, method=[envelope, range=-1..1]) ):

 

[HFloat(-3.747279823038963e-4), HFloat(0.28029046512049793)]

 

memory used=507.56KiB, alloc change=0 bytes, cpu time=7.00ms, real time=7.00ms, gc time=0ns

 

# A lengthy way

p := Probability(C < -1, numeric);
q := 1-Probability(C > 1, numeric);

U_sample := Sample(Uniform(p, q), N):

C_cdf  := unapply(CDF(C, z), z):
C_icdf := unapply(solve(C_cdf(u)=z, u), z);

C_sample := C_icdf~(U_sample):
Histogram(C_sample);

[Mean, Variance](C_sample);

CodeTools:-Usage( C_icdf~(U_sample) ):

HFloat(0.25)

 

HFloat(0.75)

 

proc (z) options operator, arrow; -cot(z*Pi) end proc

 

 

[HFloat(-0.008152438124121577), HFloat(0.2737286335656307)]

 

memory used=2.21MiB, alloc change=0 bytes, cpu time=30.00ms, real time=30.00ms, gc time=0ns

 

 


 

Download CauchySampling.mw

Completing Carl Love's reply: if you want to pick "without replacement" and "without order constraints" (that is the nth picked number is not necessarily greater than the previous one) P numbers from a list L you must use combinat:-randperm
 

restart

L := [$1..100]:

P:= 3:
combinat:-randperm(L)[1..P]

[13, 71, 67]

(1)

# compare this

for k from 1 to 10 do
  combinat:-randcomb([$1..5], 5)
end do;

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

(2)

# to this

for k from 1 to 10 do
  combinat:-randperm([$1..5], 5)
end do;

[2, 4, 5, 1, 3]

 

[5, 2, 1, 4, 3]

 

[5, 1, 2, 3, 4]

 

[5, 1, 3, 4, 2]

 

[2, 4, 1, 3, 5]

 

[2, 4, 5, 1, 3]

 

[5, 4, 2, 3, 1]

 

[1, 5, 2, 3, 4]

 

[1, 2, 3, 4, 5]

 

[1, 3, 2, 5, 4]

(3)

 


 

Download pick.mw

If you really want to use colorscheme you can use this (the key point is style=point)

N  := numelems(x1):
xy := convert(Matrix(2, N, [x1, y1])^+,listlist):

plot(xy, style=point, colorscheme=["zgradient", ["Orange", "Red", "NavyBlue"]]);  # for instance


Another way to blavk/Blue coloring:

with(plots):
display(seq(pointplot([x1[k], y1[k]], color=piecewise(z1[k]<0.5, "Black", "RoyalBlue")), k=1..N))

 


 

restart

fe := Int(sqrt(1+(-5.557990765*sin(5.557990765*x)-7.3*cos(5.557990765*x)-5.6*sinh(5.557990765*x)+7.3*cosh(5.557990765*x))^2), x = 0 .. al)

Int((1+(-5.557990765*sin(5.557990765*x)-7.3*cos(5.557990765*x)-5.6*sinh(5.557990765*x)+7.3*cosh(5.557990765*x))^2)^(1/2), x = 0 .. al)

(1)

# this to have a broad idea of a search interval

plot(evalf(subs(al=a, fe)), a=0..1, gridlines=true)

 

g := proc(a)
  evalf(subs(al=a, fe))-0.5
end proc;

proc (a) evalf(subs(al = a, fe))-.5 end proc

(2)

fsolve(g(al), al=0..1);

.1550088024

(3)

evalf(subs(al=%, fe))

.5000000002

(4)

 


 

Download int.mw

@radaar   
@Carl Love

What do you mean by "surrogate optimization built in function"?
A surrogate can be as simple as the linear model Statistics:-LinearFit returns, or as complex as a neural network (Maple 2018 and later) with, in between GPE (gaussian process emulation such as kriging for instance).
Surrogate-based optimization generally uses GPE

As Carl wrote you could use the Kriging interpolation (I think this requires Maple 2018 or later).
But this implementation of the kriging method is more related to geostatistics than to surrogate modeling of an expensive code (mainly: correlation functions not well adpated).

There also exists a (commercial) toolbox dubbed "GlobalOptimization" which  is supposed to implement the EGO (Efficient Global Optimization method based on kriging surrogates).... unfortunately it gives poor results and doesn't work as claimed (in particular it can't be parameterized as described in the help pages; I posted a question about it a few weeks ago).

My personal (and thus questionable) opinion is that Maple is not yet adpated to do Surrogate-based optimization.

I can't picture out what your problem really is...

  1. Were did you find this procedure "print", is it one you wrote somewhere else?
  2. Use Matrix instead of matrix with strng type elements.
  3. Maybe read the help pages about Matrix and print?

 


 

restart:

problems := Matrix(5, 1, [M1, M2, M3, M4, M5])

problems := Matrix(5, 1, {(1, 1) = M1, (2, 1) = M2, (3, 1) = M3, (4, 1) = M4, (5, 1) = M5})

(1)

# g := "negative effect on the health of the consumer"

PA := "probability occurrence of each problems"

"probability occurrence of each problems"

(2)

g := Matrix(1, 3, [1, 2, 3])

g := Matrix(1, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 3})

(3)

Criteria := Matrix(1,3,[ "not serious",  "moderately serious", "very serious"]);

#` For each value of g we print the corresponding adjective from the matrix Criteria `

for i from 1 by 1 to 3 do
  if  g(i)=i  then   
    print(Criteria(i));  
  end if   
end do    
 

Criteria := Matrix(1, 3, {(1, 1) = "not serious", (1, 2) = "moderately serious", (1, 3) = "very serious"})

 

"not serious"

 

"moderately serious"

 

"very serious"

(4)

IN := Matrix(1, 3, [1, 2, 3])

IN := Matrix(1, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 3})

(5)

A := Matrix(1, 3, ["without risk", "occasional dange", "very common danger"])

A := Matrix(1, 3, {(1, 1) = "without risk", (1, 2) = "occasional dange", (1, 3) = "very common danger"})

(6)

for i to 3 do
  if IN(i) = i then
    print(A(i))
  end if
end do

"without risk"

 

"occasional dange"

 

"very common danger"

(7)

 


 

Download DontUnderstandYourProblem.mw

PS: if my answer is completely out of line, please forget it 

An alternative:

For "Tally" see the corresponding help page. 
It's ouput is a list of equalities, the lhs being the "name" of an element of Set and the rhs its number of occurences.
 

Set:={{1, 2}, {1, 3}, {1, 4}, {2, 3}};

Statistics:-Tally(op~([Set[]]))

[1 = 3, 2 = 2, 3 = 2, 4 = 1]

(1)

# works also for names

Set:={{a, b, 2}, {x, 1}, {a, x}, {4, 1, a, 3}};

Statistics:-Tally(op~([Set[]]))

{{1, x}, {a, x}, {2, a, b}, {1, 3, 4, a}}

 

[1 = 2, 2 = 1, 3 = 1, 4 = 1, a = 3, x = 2, b = 1]

(2)

# and even in more general cases

Set:={{sin(a), b, 2}, {x, 1}, {a+x, x}, {4, 1, a, 3}};

Statistics:-Tally(op~([Set[]]))

{{1, x}, {x, a+x}, {2, b, sin(a)}, {1, 3, 4, a}}

 

[1 = 2, 2 = 1, 3 = 1, 4 = 1, a = 1, sin(a) = 1, x = 2, b = 1, a+x = 1]

(3)

 


 

Download tally.mw

I do not think you will be able to obtain an analytic expression.
Here is an approximate one:


 

restart:

z := sqrt(G*(2-G))+(1-G)*(arccos(G-1)+(1/5)*Pi)+k*sqrt(G*(2-G))*cos(sqrt((1+k)/k)*arccos(G-1)+(1/5)*Pi)+sqrt(k/(1+k))*(1-G)*k*sin(sqrt((1+k)/k)*arccos(G-1)+(1/5)*Pi)

(G*(2-G))^(1/2)+(1-G)*(arccos(G-1)+(1/5)*Pi)+k*(G*(2-G))^(1/2)*cos(((1+k)/k)^(1/2)*arccos(G-1)+(1/5)*Pi)+(k/(1+k))^(1/2)*(1-G)*k*sin(((1+k)/k)^(1/2)*arccos(G-1)+(1/5)*Pi)

(1)

p := plots:-contourplot(z, G=-1..2, k=1..10, contours=[0], grid=[400, 400]):

op(-1, p);
c := op(1..-2, op(1, p)):
M := convert(op~(1, [c]), Matrix);

"AXESLABELS[Typesetting:-mi("G",italic = "true",mathvariant = "italic"), Typesetting:-mi("k",italic = "true",mathvariant = "italic")]"

 

Matrix(%id = 18446744078394417014)

(2)

plots:-display(p, gridlines=true)

 

# here is the approximate relation k=f(G) such that z(k=f(G), G) = 0)

fit := Statistics:-LinearFit([1, G, G^2, G^3], M, [G]);

plots:-display(p, plot(fit, G=min(M[..,1])..max(M[..,1]), color=blue, thickness=7, transparency=0.8, gridlines=true))

HFloat(5.167125933526261)-HFloat(6.356106757782967)*G+HFloat(3.0043651079986984)*G^2-HFloat(0.767811803757935)*G^3

 

 

zg := eval(z, k=fit):  # = z(k=fit(G), G):

fit_error := evalf[4](map(u -> eval(zg, G=u), [seq(0..1, 0.01)]))

[0.15e-1, 0.3e-2, 0.5e-2, 0.2e-2, 0.2e-2, 0., 0.2e-2, -0.4e-2, -0.3e-2, -0.3e-2, -0.2e-2, 0.1e-2, -0.4e-2, -0.4e-2, -0.1e-2, -0.4e-2, -0.2e-2, 0.5e-3, -0.27e-2, -0.16e-2, -0.35e-2, -0.6e-3, 0., 0.2e-3, -0.16e-2, 0.8e-3, 0.5e-3, -0.22e-2, -0.6e-3, -0.20e-2, 0.8e-3, 0.10e-2, -0.5e-3, 0.16e-2, -0.6e-3, 0.42e-2, 0.42e-2, 0.30e-2, 0.95e-3, 0.309e-2, 0.246e-2, 0.280e-2, 0.157e-2, 0.433e-2, 0.3612e-2, 0.191e-2, 0.247e-2, 0.188e-2, 0.351e-2, 0.24e-3, 0.256e-2, 0.295e-2, 0.275e-2, 0.299e-2, 0.65e-3, 0.271e-2, 0.229e-2, 0.58e-3, 0.40e-3, 0.24e-2, 0.12e-2, 0.19e-2, 0.5e-3, 0.6e-3, 0.10e-2, -0.2e-3, -0.6e-3, -0.10e-3, -0.92e-3, -0.98e-3, -0.19e-3, -0.199e-2, -0.288e-2, -0.349e-2, -0.258e-2, -0.166e-2, -0.225e-2, -0.358e-2, -0.265e-2, -0.176e-2, -0.248e-2, -0.318e-2, -0.339e-2, -0.304e-2, -0.336e-2, -0.198e-2, -0.170e-2, -0.224e-2, -0.183e-2, -0.236e-2, -0.79e-3, -0.107e-2, -0.133e-2, -0.57e-3, 0.30e-3, 0.29e-3, 0.140e-2, 0.1598e-2, 0.2962e-2, 0.3399e-2, 0.4400e-2]

(3)

 


 

Download LevelCurve2.mw

Here is a solution where each level curve is a colred differently

Maybe you would also have a legend for each level curve?
In this case you could ne interested by what acer did here 2D contour plot and legend

PS I removed the ouput for the uploaded display doesn't have colored level curves (but they are)
 

restart:

with(plots):

f := sin(x*y):
r := x = -2 .. 2, y = -2 .. 2:

cols := [seq(ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]),i=1..19)]:
plots:-display(
  plot3d(f, r, axes = framed, style=surface, color=gray),
  seq(contourplot3d(f, r, contours=[(j-10)/10], color=cols[j]), j=1..19)
);  

 

 


 

Download contourplot.mw

Luckily I had already done it, I just wasted a little time finding the code.
From memory I believe I recoded in Maple the code from R (not completely sure).

Here it is. 
I remember having done some tests and it seemed ok by the time , but feel free to inform me if it's not.


 

restart:

LSD_Sampling := proc(alpha, N)
  local c, V, X, n, u, q:

  uses Statistics:

  c     := log(1-alpha):
  V     := Sample(Uniform(0, 1), N):
  X     := Vector[column](N):
  for n from 1 to N do
     if V[n] > alpha then
        X[n] := 1:
     else
        u := Sample(Uniform(0, 1), 1)[1]:
        q := 1 - exp(c*u):
        if V[n] < q^2 then
           X[n] := ceil(1+log(V[n])/log(q)):
        elif V[n] < q then
           X[n] := 1:
        else
           X[n] := 2:
        end if:
     end if:
  end do:

  return X:
end proc:

Statistics:-Histogram(LSD_Sampling(0.9, 10000))

 

 


 

Download LogarithmicSeriesSampling.mw


BTW: I'm curious, what application of the Logarithmic Series distribution are you interested in?
 

Even Re(code) doesn't work. The reason is: Maple doesn't know if the variables are real or complex !

Assuming all the quanties are real a simple way is

subs(J=I, collect(subs(I=J, code), J));

Other possibility: use assume, or assuming (see help pages)

I already ask this question a few years ago, you can trace the answers here

Question: I'd like to have some tips about Version Control (source control) process

The solution I use:

  1. I assume you have 2 mw files, let's say File1.mw and File2.mw whoch correspond to different versions of the same code.
     
  2. Open File1.mw and Export it to maple input format (I choose the same name as the mw file)
     
  3. Open File2.mw and Export it to maple input format; two files File1.mpl and File2.mpl are created
     
  4. If not already installed, download NotePad++ (free on Windows ; Notepad++ has contextual coloring for Maple syntax)
     
  5. Open NotePad++ and load File1.mpl and File2.mpl
     
  6. In the menu of NotePad++ choose Run > Compare (this from memory because I do not have NotePad++ on this laptop)
    The differences between the two versions are then easily visble

Once done, you have the choice to keep developping with the Java interface of Maple (as usual), or to develop "within" NotePad++.
In this latter case you will have to load an mpl file in a Maple session to verify that its content gives the expected result.
It's up to you, I personally prefer developping "within" a worksheet, saving in a new mw file, and exporting it in a new mpl file.

PS : there exist probably tools apart NotePad++ that can "understand" Maple's syntax and propose contextual coloring.

I would like to complete acer's answer. Personally I like to use notations like theta__i instead of theta[i], because i is not only limited to integer values.
But using theta__i instead of theta[i] requires some precautions.
You will find below a few examples.

Be careful, when used in a procedure, variables of the form theta__||n (last example) are not implicitely considered as local but as global: this may cause unwanted behaviours.

 

# differences :

a := theta[2]:
b := theta__2:

whattype(a), whattype(b);

printf("%a, %a", theta[2], theta__2);
print(a,b);

indexed, symbol

 

theta[2], theta__2

 

theta[2], theta__2

(1)

# use of theta__  ; first case: wrong way

restart:

vars := seq(theta__i, i=1..4);

mean := Array(1..4, i -> diff(w(vars), theta__i));

 

vars := `#msub(mi("&theta;",fontstyle = "normal"),mi("i"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("i"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("i"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("i"))`

 

Array(%id = 18446744078232662014)

(2)

# use of theta__  ; second case: good way

restart:

vars := seq(theta__||i, i=1..4);

mean := Array(1..4, i -> diff(w(vars), theta__||i));

 

vars := `#msub(mi("&theta;",fontstyle = "normal"),mi("1"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("2"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("3"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("4"))`

 

Array(%id = 18446744078232662014)

(3)

# Why using theta__ instead of theta[..] can be useful?
# Because the index n in theta__n can be of non numeric type:

restart:
MyIndices := [a, b, c, d];
vars := seq(theta__||n, n in MyIndices);

mean := Array(1..4, i -> diff(w(vars), vars[i]));

MyIndices := [a, b, c, d]

 

vars := `#msub(mi("&theta;",fontstyle = "normal"),mi("a"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("b"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("c"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("d"))`

 

Array(%id = 18446744078232662014)

(4)

restart:

f := proc()
local n:
for n from 1 to 4 do
  theta__||n := n
end do:
end proc:

f():
seq(theta__||i, i=1..4);

1, 2, 3, 4

(5)

 


 

Download indexed_theta.mw

 

Two ways (probably among others)

restart:
vars := seq(theta[i], i=1..4);
mean := Array(1..4, i -> diff(w(vars), theta[i]));

# or
mean := Array(1..4, i -> diff(w(vars), vars[i]));

Download mean.mw

In case you would be interested, give a lookk to this site :

https://oeis.org/?language=english

Next just copy-paste 2,3,4,6,8,9,10,12,14,15,16,18 in the search field and enjoy

 

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