mmcdara

7891 Reputation

22 Badges

9 years, 59 days

MaplePrimes Activity


These are replies submitted by mmcdara

@acer  @vv

Thanks for the answer.

What I did seems simpler (maybe it has some drawbacks when compared to your solution?)
X := (mu, sigma) -> RandomVariable(Normal(mu, sigma)):
Mean(X(0, 1));
PDF(X(0, 1), t)
;

 

I've just got it : its name is "Specialize", from Statistics package !!!

Thanks to both of you for your involvment

@vv  
Pretty workaround.
But you need to write Mean(eval(X)) to get mu.

Too bad I can't find the command I saw a few times ago.

@Carl Love 

Right.
acer eventually gave me the solution, but I nevertheless sent a reply to dharr to explain the motivation of my question.

 

@dharr 

You're right, given the pdf or the cdf one can recognize the distribution law.
But what I was looking for was different (acer gave me the answer and this is why I use a oaste tense).

To give you idea of my motivation, please look to the Statistics[ChiSquareGoodnessOfFitTest] help page.
In the example (Maple 2015.2) the "empirical distribution" Ob is compared to an "exact distribution" Ex. In this example Ex is obviously a representation of some uniform (discrete or not) random variable.
Suppose now that you have a sample S drawn from some population P of unknown distrinution D and that you infer that the distribution of P is D*.
To check the null hypothesis "S is drawn from a population distributed according to D* " you can use this Chis square goodness of fit test. But as it's written it operates on two vectors containing the numbers of occurrences in different classes.

What I want to do is to build a procedure Chi2GoFTest with 2 parameters S and D* which returns the probability to reject wrongly this null hypothesis.

Something like 
PriorDistribution := Statistics:-RandomVariable( blablabla );
Chi2GoFTest(S, PriorDistribution);

Nothing complicated here.
But I would like this procedure prints lines of the form:
Null hypothesis: "the sample S is drawn from a distribution blablabla"
Probability to reject wrongly the null hypothesis = xxxx


When I call Chi2GoFTest the "object" PriorDistribution is passed and I can manipulate it within this procedure, for instance I can print its pdf or compute statistics or generate samples. But I didn't know how to print blablabla.
Eventually acer gave me the answer.

Thanks for your  involvement

@acer 

Thanks acer, it's exactly what I was looking for.

@radaar 

Give the terms you use it seems to me you are familiar with the topic.
I guess your question could be summarize this way:

  • I have a "code" (for instance the procedure dsolve (..., numeric) can return) wich relates some inputs a single (to simplify) quantity of interest.
    I would like to find an input vector such that this quantity of interest (qoi) reaches its minimum value.
    Unfortunately my "code" is rather expensive to evaluate (so I can't do intensive blind search to assess this input vector), and I don't know closed form expressions of the derivatives of the qoi relative to each input (which avoids using some minimization algorithms).
     
  • So the idea to realize a moderate number of runs of the code for an a priori set of inputs; to infer an approximate model  which would represent the variations of the qoi given the inputs; at last to use this model as a surrogate of the original "code".


I may be completely out of line but, if not, the first problem relates to the reliability of the surrogate: how well does it represent the "code". Obviously a poor surrogate will give you a poor minimizer of the qoi.

This surrogate can be either local or global.  
Kriging is aimed to build global surrogates and there exist a widely used strategy (Expectation Improvment) used to enhanced iteratively the quality of the Kriging surrogate (Expectation improvement is a global optimization problem where the qoi represents the improvement of the representativity of the surrogate). Once a reliable surrogate has been obtained you can use it in a global optimization method to minimize your original qoi.
If you are interested in kriging based optimization the most efficient procedure is dubbed Efficient Global Optimization procedure (EGO).

Surrogate based optimization using other types of surrogates (neural networks, regression models and so on) are generally more unstable and reliable.

Using local (linear or quadratic) surrogates is also a way you could use.
I could advice you this book: Derivative-free and blackbox optimization, Charles Audet, Warren Hate (Springer) and in particular Part 4: Model based methods.

But in any way you will have to write a lot of code by yourself (the book clearly describes the algorithms).

Good luck

 

@MDD 

Perhaps this question about afficiency is of major concern to some people on this site. But I'm not one of them. So I'll be happy if you prefer my answer... but not at all sad if you'd prefer Acer's.

@acer 

Did you change your answer meanwhile?

If not I'm sorry, this means I stopped reading your answer just after the first group of commands.... 

 

@acer 

I'm not sure it is possible to answer the question correctly without explicitely declaring what the indeterminates are.
Following the OP you explicitely considered they were x and y.
But what makes the difference between F[1] and F[2] when F := [x, a]?

I believe a more reliable way is to declare the indeterminates first.
Then, assuming each member of F is a polynomial, a simple way to determine if F[n] is a constant polynomial relatively to the declared indeterminates, is to check if its degree is equal to 0.

 

restart:

F := [a*x^3-1, b*y, -1, c, b+a];
V := [x, y]:
map(degree, F, V);

[a*x^3-1, b*y, -1, c, b+a]

 

[3, 1, 0, 0, 0]

(1)

F := [a*x^3-1, b*y, -1, c, b+a];
V := [a, c]:
map(degree, F, V);

[a*x^3-1, b*y, -1, c, b+a]

 

[1, 0, 0, 1, 1]

(2)

 


 

Download map-degree.mw

@Carl Love 

Hi, I'm not sure your answer is complete (maybe I didn't understand correctly what the question really was).
It seems to me that the OP wants to apply its "removing rule" iteratively?

For instance


L := $1..100]:
M := RemoveEveryNth(L, 2);

k := 1:
while M <> L do
  k := k+1:
  R := ({M[]} intersect {L[]})[k] ;
  L := copy(M):
  M := RemoveEveryNth(L, R):
  print(M);
end do:

@BjaMaple 

Also : seq(x__||i, i=1..4);

@Rouben Rostamian  

Given that  Xbar = Sum(x[i], i=1..n)/n one easily finds that diff(Xbar, x[i]) = 1/n whatever the value of i
Then diff(b, x[i]) = diff(b, Xbar)*diff(Xbar, x[i]) = diff(b, Xbar)/n
(which is just "good sense" if we refer to the statistical meaning of Xbar).

Your procedure c_2 computes diff(b, x[i]). 
You will find that  sum(c_2(my_x, my_xbar, i), i=1..numelems(my_x)) = 0  ***(see below)
But, given diff(b, x[i]) = diff(b, Xbar)/n, applying the summation on the lhs and rhs of this equality should imply diff(b, Xbar)=0
You find a value (c_1) of -1.479...

This comes from the fact that, contrary to the requirement, x[i] and Xbar are not independent.
This hypothesis farazhedayati asked you to use is inconsistent

About the result sum(c_2(my_x, my_xbar, i), i=1..numelems(my_x)) = 0 
This result holds for any list my_x
(try for instance my_x := convert(Statistics:-Sample(Uniform(0, 1), 17), list) ).
This is a direct consequence of the fact that Xbar = Sum(x[i], i=1..n)/n

@farazhedayati 

"Please consider Xbar and x[i] as independent variables" : Wrong, as xbar=sum(x[i], i=1..n)/n they are obviously not independent
diff( sum(x[i], i=1..n)/n, x[i]) = 1/n whatever the value of i
 

@9009134 

Sorry for this late reply.

Here is an attempt to answer your question.
Be carefull: the method I used in my previous post can give false results if the integrand is singukar over the integration domain (look yellow highlighted red text)


 

restart

# some examples

  a := diff(f(x,y), x$3) + diff(f(x,y), [x, y]) + diff(f(x,y), y$2) + f(x, y) + diff(g(x, y), y) + h(y);
# a := diff(f(x,y), x$3) + diff(f(x,y), [x, y]);
# a := f(x, y) + h(y);

diff(diff(diff(f(x, y), x), x), x)+diff(diff(f(x, y), x), y)+diff(diff(f(x, y), y), y)+f(x, y)+diff(g(x, y), y)+h(y)

(1)

# define the alias ExI

alias(ExI=expand@int);

# hasfun(a, diff) returns true if expression a contains terms with derivatives

if hasfun(a, diff) then
  b  := [op(select(has, a, diff))];                 # extract all terms with derivatives

  c  := op~(-1, [op(%)]);                           # extract the last differentiation variable for each term in b
                                                    # c will be used to integrate with respect to these last
                                                    # differentiation variables
                   
  d  := map(C -> if C=x then y else x end if, c);   # d denotes the other integration variable

  zip((u, v) -> int(u, v), b, c);                   # integrate ech term in b along its last differentiation variable
  I1 := add(zip((u, v) -> int(u, v), %, d));        # integrate each previous term along the second integration variable
end if:

# add to I1 the double integrations with respect to x and y of all the terms that do not contain derivatives

I2 := I1 + ExI(ExI(select(`not`(has), a, diff), x), y)

ExI

 

int(diff(diff(f(x, y), x), x), y)+f(x, y)+int(diff(f(x, y), y), x)+int(g(x, y), x)+int(int(f(x, y), x), y)+x*(int(h(y), y))

(2)

# boundaries of the integration domain Omega (cartesian prodict of BX by BY)

BX := [x__1, x__2];
BY := [y__1, y__2];

[x__1, x__2]

 

[y__1, y__2]

(3)

# example

I3 := eval(I2, [f(x, y)=x*cos(y), g(x, y)=exp(x*y), h(y)=1/y^3])

x*cos(y)+exp(x*y)/y-(1/2)*x/y^2

(4)

# matrix of the evaluations of I3 at the 4 corners of Omega

M := Matrix(2, 2, (i, j) -> eval(I3, [x=BX[i], y=BY[j]]))

M := Matrix(2, 2, {(1, 1) = `#msub(mi("x"),mi("1"))`*cos(`#msub(mi("y"),mi("1"))`)+exp(`#msub(mi("x"),mi("1"))`*`#msub(mi("y"),mi("1"))`)/`#msub(mi("y"),mi("1"))`-(1/2)*`#msub(mi("x"),mi("1"))`/`#msub(mi("y"),mi("1"))`^2, (1, 2) = `#msub(mi("x"),mi("1"))`*cos(`#msub(mi("y"),mi("2"))`)+exp(`#msub(mi("x"),mi("1"))`*`#msub(mi("y"),mi("2"))`)/`#msub(mi("y"),mi("2"))`-(1/2)*`#msub(mi("x"),mi("1"))`/`#msub(mi("y"),mi("2"))`^2, (2, 1) = `#msub(mi("x"),mi("2"))`*cos(`#msub(mi("y"),mi("1"))`)+exp(`#msub(mi("x"),mi("2"))`*`#msub(mi("y"),mi("1"))`)/`#msub(mi("y"),mi("1"))`-(1/2)*`#msub(mi("x"),mi("2"))`/`#msub(mi("y"),mi("1"))`^2, (2, 2) = `#msub(mi("x"),mi("2"))`*cos(`#msub(mi("y"),mi("2"))`)+exp(`#msub(mi("x"),mi("2"))`*`#msub(mi("y"),mi("2"))`)/`#msub(mi("y"),mi("2"))`-(1/2)*`#msub(mi("x"),mi("2"))`/`#msub(mi("y"),mi("2"))`^2})

(5)

# definite integration of I3 over Omega

I4 := M[2, 2] + M[1, 1] - M[1, 2] - M[2, 1]

x__2*cos(y__2)+exp(x__2*y__2)/y__2-(1/2)*x__2/y__2^2+x__1*cos(y__1)+exp(x__1*y__1)/y__1-(1/2)*x__1/y__1^2-x__1*cos(y__2)-exp(x__1*y__2)/y__2+(1/2)*x__1/y__2^2-x__2*cos(y__1)-exp(x__2*y__1)/y__1+(1/2)*x__2/y__1^2

(6)

# verification

A := eval(a, [f(x, y)=x*cos(y), g(x, y)=exp(x*y), h(y)=1/y^3]);

# WATCHOUT : the above strategy doesn't handle possible singularities inside Omega
#
# ILLUSTRATION
ia := int(int(A, x=BX[1]..BX[2]), y=BY[1]..BY[2]);

# I cheat here

IA := int(int(A, x=BX[1]..BX[2]), y=BY[1]..BY[2]) assuming BX[1] < BX[2], 0 < BY[1], BY[1] < BY[2]

-sin(y)+x*exp(x*y)+1/y^3

 

Warning, unable to determine if 0 is between y__1 and y__2; try to use assumptions or use the AllSolutions option

 

int((x__1*sin(y)*y^3-x__2*sin(y)*y^3-exp(x__1*y)*x__1*y^2+exp(x__2*y)*x__2*y^2+exp(x__1*y)*y-exp(x__2*y)*y-x__1+x__2)/y^3, y = y__1 .. y__2)

 

-(1/2)*(2*cos(y__2)*x__1*y__1^2*y__2^2-2*cos(y__2)*x__2*y__1^2*y__2^2-2*cos(y__1)*x__1*y__1^2*y__2^2+2*cos(y__1)*x__2*y__1^2*y__2^2+2*exp(x__1*y__2)*y__1^2*y__2-2*exp(x__2*y__2)*y__1^2*y__2-2*exp(x__1*y__1)*y__1*y__2^2+2*exp(x__2*y__1)*y__1*y__2^2-x__1*y__1^2+x__1*y__2^2+x__2*y__1^2-x__2*y__2^2)/(y__1^2*y__2^2)

(7)

simplify(IA - I4)

0

(8)

# BE CAREFULL ....

int(diff(1/x, x), x=-1..1);

# but

int(diff(F(x), x), x);
eval(%, F(x)=1/x);
eval(%, x=1) - eval(%, x=-1)

-infinity

 

F(x)

 

1/x

 

2

(9)

 


 

Download Sketch.mw



 

@vv 

You're right
 

restart:

L := (r, theta) -> -2* diff(f(r, theta), r$2)+12*f(r, theta)-12*diff(f(r, theta), theta$2)/r+0.6*(diff(f(r, theta), theta, theta, r, r))*(1/4)+0.5*(diff(f(r, theta), theta, theta, theta, theta))*(1/4)

proc (r, theta) options operator, arrow; -2*(diff(f(r, theta), `$`(r, 2)))+12*f(r, theta)-12*(diff(f(r, theta), `$`(theta, 2)))/r+.6*(diff(f(r, theta), theta, theta, r, r))*(1/4)+.5*(diff(f(r, theta), theta, theta, theta, theta))*(1/4) end proc

(1)

S := L(r, theta);

-2*(diff(diff(f(r, theta), r), r))+12*f(r, theta)-12*(diff(diff(f(r, theta), theta), theta))/r+.1500000000*(diff(diff(diff(diff(f(r, theta), r), r), theta), theta))+.1250000000*(diff(diff(diff(diff(f(r, theta), theta), theta), theta), theta))

(2)

IS := int(int(op(1, S), r), theta) + int(int(op(2, S), r), theta) + int(int(op(3, S), theta), r) + int(int(op(4, S), theta), r) + int(int(op(5, S), theta), r)

int(-2*(diff(f(r, theta), r)), theta)+int(int(12*f(r, theta), r), theta)+int(-12*(diff(f(r, theta), theta))/r, r)+.1500000000*(diff(diff(f(r, theta), r), theta))+int(.1250000000*(diff(diff(diff(f(r, theta), theta), theta), theta)), r)

(3)

F := (sin(-4.700000000*10^(-6)+4.700000000*r)-0.1369508410*sinh(-4.700000000*10^(-6)+4.700000000*r))*cos(6*theta);

(sin(-0.4700000000e-5+4.700000000*r)-.1369508410*sinh(-0.4700000000e-5+4.700000000*r))*cos(6*theta)

(4)

subs(f(r, theta)=F, IS):

# Exact primitive

ISF := eval(%);

.1666666667*(-9.400000000*cos(-0.4700000000e-5+4.700000000*r)+1.287337905*cosh(-0.4700000000e-5+4.700000000*r))*sin(6.*theta)+.1666666667*(-2.553191489*cos(-0.4700000000e-5+4.700000000*r)-.3496617217*cosh(-0.4700000000e-5+4.700000000*r))*sin(6.*theta)+(0.3384000000e-3+0.*I)*sin(6.*theta)*Ei(1., -(4.700000000*I)*r)+4.930207104*sin(6.*theta)*Ei(1., -4.700000000*r)-4.930253448*sin(6.*theta)*Ei(1., 4.700000000*r)+(-113.0973355-0.5315574770e-3*I)*sin(6.*theta)*csgn(r)+(72.00000000+0.3384000000e-3*I)*sin(6.*theta)*Si(4.700000000*r)-.9000000000*(4.700000000*cos(-0.4700000000e-5+4.700000000*r)-.6436689527*cosh(-0.4700000000e-5+4.700000000*r))*sin(6*theta)+0.2700000000e-7*sin(6.*theta)*(-212765957.4*cos(-0.4700000000e-5+4.700000000*r)-29138476.81*cosh(-0.4700000000e-5+4.700000000*r))

(5)

Ar := [0, 1]:
At := [0, 2*Pi]:

sISF := simplify(ISF):
M := Matrix(2, 2, (i,j) -> evalf(eval(limit(sISF, r=Ar[i]), theta=At[j])));

M := Matrix(2, 2, {(1, 1) = Float(undefined), (1, 2) = Float(undefined), (2, 1) = 0.+0.*I, (2, 2) = 0.+0.*I})

(6)

IL := M[1, 1] + M[2, 2] - M[1, 2] - M[2, 1]


 

Download UndefinedIntegral.mw

 

 

First 123 124 125 126 127 128 129 Last Page 125 of 154