mmcdara

640 Reputation

10 Badges

3 years, 44 days

MaplePrimes Activity


These are answers submitted by mmcdara

I'm not sure I understand you perfectly.
What I think is: you have a discrete RV "X" with mass function fX(x)=1/(x+1) for x = 0, ..., k, where k is some given number (>0).

If It is that, a way to proceed is to use a ProbabilityTable Distribution

 

restart:

with(Statistics):

# Discrete distribution from a probability table
# The sum of all masses must be equal to 1.

P := k -> [seq(1/(x+1), x=0..k)] /~ add(1/(x+1), x=0..k) ;

# example:

K := 5:
N := 10^3:

P(K);

X := RandomVariable(ProbabilityTable(P(K))):

S := Sample(X, N):
S[1..10];
Histogram(S);

 

P := proc (k) options operator, arrow; `~`[:-`/`]([seq(1/(x+1), x = 0 .. k)], ` $`, add(1/(x+1), x = 0 .. k)) end proc

 

"[[Typesetting:-mfrac(Typesetting:-mn("20",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("10",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("20",mathvariant = "normal"),Typesetting:-mn("147",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("5",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("4",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("10",mathvariant = "normal"),Typesetting:-mn("147",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false")]]"

 

Vector[row]([3., 2., 6., 2., 1., 1., 2., 1., 2., 1.])

 

 

 


 

Download DiscreteRV.mw

Hi, 

If your only concern is to compare the model and the experiment, and if this experiment is a "not too large" collection of triples (t, x, y) , then the option output=Array[ ...], where ... is the list of the experimental t values, is probably the simpler solution.

Regarding now your "comparison metric".
The answer is very simple: there is unfortunately no good way to proceed!
Computing "the Root Mean Square of all (xexperiment(t)-xmodel(t), yexperiment(t)-ymodel(t))" seems very reasonable.
But it can be a very poor solution for some problems, specifically if x or y are of different natures, vary on very different scales and so on.
Luckily it doesn't seem to be the case here.

Another point you must be aware of: if the (t, x, y) experimental values come from measurements, there are probably entailed by "errors". Suppose that x and y are perfectly measured, but that t is not.
For instance the "true" t is T=1, and the measured value is t=1.1: then comparing the experiment at time t=1.1 to the model at time t=1.1 doesn't make sense.
A (the?) correct approach in this case is to use a bayesian framework where one infer the true value of T given its measure t.

Maybe you could be interested in "fitting" your model on the experiments? For instance to assess the value of A ...?
You also mention that "The problem is however that even if the (x,y) curves are rather close, if may not be a good model".
Are you concerned with aapossible "model error"?
Here again, for these two questions, the correct answer relies to the "bayesian calibration of computer model" issue.

Let me know if you are interested by taht.

To conclude: try your  "Root Mean Square of all (xexperiment(t)-xmodel(t), yexperiment(t)-ymodel(t)", examine carefully the residues (the model-to-experiment distance at each time t) and if all looks good, then it's done 
 

MAT := < header1, header2, <Column(t1[2,1],1..2) | Column(t2[2,1],2) > >:
interface(rtablesize=max(Dimensions(MAT))):
MAT;

 

 

Assuming the DEFINITION below suits you, you will find some elements in the attached file.

DEFINITION:

  • Let E the ellipse defined by the equation (x/a)^2+(y/b)^2=1
  • For each point P on E one define s(P) as the measure of the angle (counter clockwise counted) between Ox and OM
  • Let A and B two points on E such that
    • s(A) > =  0
    • s(B) > s(A)
    • s(B) <= 2*Pi
  • Let M any point such that s(A) <= s(M) <= s(B)
  • Let x(A) the abscissa of A and x(B) the one of B

Then the mean value of the distance d(O, M) from O=(0, 0) to M as M moves from A to B is defined by the integral of d(O, M) between x(A) and x(B), divided by x(B)-x(A).

In file Ellipse2.mw you will find a procedure 'md' which computes the mean value of d(O, M). It takes 4 arguments:

  1. b
  2. alpha = s(A)
  3. beta = s(B)

'md' handles (if I'm not mistaken) the 3 cases:

  1. alpha < Pi and beta < Pi
  2. alpha < Pi and beta > Pi
  3. alpha > Pi and beta > Pi 


Before that you will find a few lines that treat more formally the first situation alpha < Pi and beta < Pi

As claimed in this file, I was not capable to obtain a general formal solution for any values of a and b (that's why I arbitrarily wrote a:=2 and b:=1 at the top at the file).
The search of a general formal solution requires using assumptions about a, b, p, q (see the commands in the file) and the 'int' procedure of Maple is not very comfortable with 'assume' or 'assuming' statements (someone here will probably correct me on this point)

Ellipse2.mw

Here is a possible (NOT THE) definition of the mean radius.

I replace "arithmetic mean" by mean in the sense (one sense) mentioned by Rouben:

  • I consider the canonical ellipse with principal axis Ox and Oy (natural cartesian coordinates)
  • let M1, M2<> M1 two points  on the ellipse
  • Let M a third point on the ellipse.
     Assuming M moves from M1 to M2 counterclockwise in such a way that the angular sector swept by M is less than 2*Pi,
    one can define the mean radius by the mean value of OM when M moves from M1 to M2.

Ellipse.mw

 

 

Maybe the shortest solution?
(just for fun, no practical interest when compared to the previous answers)

restart:

splitChange := proc(w)
  local bw:
  bw := map(u -> if w[u+1]<>w[u] then cat(w[u], " ") else w[u] end if, [$1..length(w)-1]);
  StringTools:-Words(cat(op(bw), w[-1]));
end proc:
  
splitChange("WPKCPYWYYYXHYY")


Download Shorter.mw

I guess you took as a reference the article by Faith Chelimo Kosgei?
(https://www.journalijdr.com/sites/default/files/issue-pdf/11819.pdf  practically the one I was able to find on the subject)

If it so I understand you asked for some help (I 've rarely read such a muddled, obscure and confusing paper!)

Nevertheless, after a lot of work spent to read between the lines of the article, you will find in the attached file the algorithm applied to the 3 toy examples this paper contains.

NewtonLagrange.mw

I leave you the pleasure to apply the method to your own test cases

I don't know what Maple version you use (in a previous thread you said tou used Maple 13?) but in Maple 2015 (and beyond.. but I can verify this right now) your code is not correct.
The line

should be replaced by 
V:= r -> piecewise(0<=r and r<=R, -VV, r>R, 0):

Here is your original code

A few points beyond the above correction

  1. Some lines are unnecessary:

     
  2. As Carl Love wrote, you need another IC/BC condition, dot instance D(u)(0) = something
    With, for instance, SYS2:={diff(u(r), r, r)+k^2*u(r) = V(r)*u(r), u(0) = 0, D(u)(0)=1} you should obtain this curve
    odeplot(sol, [r, u(r)], r = 0 .. 2);

The GB is reduced w.r.t the plex order (proof here)

GB.mw

 

It seems indeed that the help pages do not mention if th GB is reduced or not.
 

 

 

Use "randomize" to generate a random sed on each node

restart:
with(Grid):
with(Statistics):

Samplor := proc(N) 
  local Seed:
  Seed:=randomize();
  Sample(Uniform(0, 1), N); 
end proc:


Set(Samplor);
Map(Samplor, [2, 2, 2, 2]);


 

 


restart:
# example
N:=10:
Y := Array(1..N, i -> i^2);
             [1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
#
DY := map(u -> Y[u+1]-Y[u-1], [$2..N-1])
                [8, 12, 16, 20, 24, 28, 32, 36]

 

Hi, 

You will find some hints in the attached file.

There are basically 2 problems wiith your code:

  1. value is a reserved word in Maple
  2. You don't read teploty.txt, so the test fails because the "value" is not valuated.
    ReadWrite.mw

A temporary solution until someone here will provide a clever one

restart:
M := `<|>`(`<,>`(1, 3, 2), `<,>`(4, 1, 3)):

# Structure of M

S := op(2, M);

 {(1, 1) = 1, (1, 2) = 4, (2, 1) = 3, (2, 2) = 1, (3, 1) = 2,  (3, 2) = 3}
 

# Reverse S

T := map(u -> rhs(u)=[lhs(u)], S);

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

# A simple procedure to find a value "n"
# (should be improved to handle the case "n is not an element of M")

Find := proc(n) map(u -> if lhs(u)=n then rhs(u) end if, T); end proc:

Find(1);
Find(2);
Find(3);
                        {[1, 1], [2, 2]}
                            {[3, 1]}
                        {[2, 1], [3, 2]}

Find.mw

 

fsolve(f, omega=4)
               4.450852611 - 1.509500329*10^(-25)*I

fsolve(f, omega=5)
               4.450852611 - 5.580307451*10^(-24)*I


As you observe the solution is complex (even the smallness of the imaginary part can suggest it is a numerical effect).
But this info is usefull and give you a hint to search a solution in a given interval:

fsolve(f, omega=1, complex, {omega=1-I..5+I}  );  # search in the xomplex square [1, 5]x[-I, I]





By the way: Do not use simplify in fsolve to spare computational time.

 

 

Your integral is just the mathematical expectation of the random variable 0.3*p1+0.7*p2:

Mean(0.3*p1+0.7*p2)
   0.01669578723

No need to look further

More than this :
CodeTools:-Usage(Mean(0.3*p1+0.7*p2))
memory used=101.54KiB, alloc change=0 bytes, cpu time=2.00ms, real time=19.00ms, gc time=0ns
                         0.01669578723


... and finally:
Mean(3/10*p1+(1-3/10)*p2)

3 4 5 6 7 8 9 Page 5 of 9