sand15

822 Reputation

13 Badges

10 years, 246 days

MaplePrimes Activity


These are answers submitted by sand15


Convolution is useful to build the PDF of the sum of two (or more) independent random variables but is of no help to compute the PDF of a function ( |.| for instance) of a random variable.

Furthermore, you write "..., but the "rules" of convolution, particularly setting the upper and lower boundaries of the integrals confuse me.".
This is very simple: as you are dealing with PDF of continuous random variables (so functions defined over R) the lower bound is minus infinity and the upper bound is plus infinity, that's all.
You can of course be more astute while integrating over the support of one of the two random variables, but I do not recommend you doing so: read Convolution_bounds.mw to see what mistakes you could do.
In fact, integrating over the smallest necessary range presents an interest only in numerical integration, but if you integrate exactly this trick has absolutely no utility.



Here is the solution in three steps.

restart


I tried to reduced to the minimum the use of Statistics built-in functions

with(Statistics):


U and V iid Uniform(0, 1)

Step 1: PDF(U-V)

Let us write nV = -V

__U  := RandomVariable(Uniform( 0, 1)):
__nV := RandomVariable(Uniform(-1, 0)):

phi[U]  := unapply(PDF(__U, z), z):
phi[nV] := unapply(PDF(__nV, z), z):
 

phi[U-V] := unapply( int(phi['U'](u)*phi['nV'](z-u), u=-infinity..+infinity), z)

proc (z) options operator, arrow; piecewise(z <= -1, 0, z <= 0, 1+z, z <= 1, 1-z, 1 < z, 0) end proc

(1)


Step 2: PDF(|U-V|)

# PDF(|U-V|) is easily obtained using CDF(U-V).
# You can use the built-in Probability function, or do elementary geometry.
#
# For instance: Prob(|U-V| < 0.7) = gray area = 1 - green areas

d := 0.7;
plots:-display(
  plot(phi[U-V](z), z=-1..1, color=blue, thickness=3)
  , plot(phi[U-V](z), z=-d..d, color=gray, filled=true)
  , plot(phi[U-V](z), z=-1..-d, color="Chartreuse", filled=true)
  , plot(phi[U-V](z), z=d..1, color="Chartreuse", filled=true)
  , scaling=constrained
)

.7

 

 

# The green area is

d := 'd':
green_area := 2 * ``((1-d)^2/2)

2*``((1/2)*(1-d)^2)

(2)

# So Prob(|U-V| < d) equals

`Prob(|U-V| < d)` = 1 - expand(green_area);

`Prob(|U-V| < d)` = 1-(1-d)^2

(3)

# Meaning the restriction of the PDF of |U-V| to 0 <= d <= 1 is

phi[`|U-V|`] := diff(rhs(%), d)

2-2*d

(4)

# The PDF of |U-V| then is

phi[`|U-V|`] := unapply(piecewise(d < 0, 0, d < 1, phi[`|U-V|`], 0), d)

proc (d) options operator, arrow; piecewise(d < 0, 0, d < 1, 2-2*d, 0) end proc

(5)


Let X and Y iid Uniform(0, 1)

Step 3: PDF( |U-V| + |X-Y| )

By convolution:

phi[`|U-V| + |X-Y|`] := unapply( simplify(int(phi[`|U-V|`](u)*phi[`|U-V|`](d-u), u=-infinity..+infinity)), d):

phi[`|U-V| + |X-Y|`](d)

piecewise(d <= 0, 0, d <= 1, (2/3)*d^3-4*d^2+4*d, d <= 2, -(2/3)*d^3+4*d^2+16/3-8*d, 2 < d, 0)

(6)

 

 

Download DensityFunction_sand15.mw


BTW: in Statistics the distance you use is named Manhattan distance or "Taxicab distance".

 


My explanation intuition

On today machines 1 hyperthreaded proc corresponds to 1 phyisical proc + 1 logical proc.
In Grid language this means 1 core = 1 physical node + 1 logical node, so setting 4 nodes means 2 cores (=2 proc)

BTW: there is a course about Grid package somewhere here (a serie of two posts I wasn't capable to put the finger on,  but I'm sure I read them years ago...)
 

@brian bovril 

I've Googled the Round-Robin method (and discovered it was what is named the  "méthode du Tourniquet" in French, see picture below)

(this is my professional login, aka mmcdara at home)

Nevertheless I didn't find anything a Frenchie can easily understand.
So here is the last state of my work

RoundRobin.mw

Result:

L := [3, 4, x, x^2, x^3, sin(x)]:
select(t -> t::numeric and t > 1, degree~(L, x))
                             [2, 3]

# All the exponents but 0 and 1
L := [3, 4, x, x^2, x^3, sin(x), x^(-4)]:
select(t -> t::numeric and not(member(t,{0, 1})), degree~(L, x))
                          [2, 3, -4]


This piece of code "extracts" all integer exponents of x but 0 and 1

L := {3, 4, x, x^2, x^3, a(x^5), x^(-4), b(a(x^6)), exp(x), b(x^7, x^8)}:

# F := select(t -> hastype(t, `^` ), L); # is redundant
F := copy(L) # is enough because op discards 3, 4, x

while indets(F minus {x}) <> {} do
  F := op~(F);
end do:
F minus {x}
                   {-4, 2, 3, 4, 5, 6, 7, 8}


In case you would like to keep exponent 1 too:

L := {3, 4, x, x^2, x^3, a(x^5), x^(-4), b(a(x^6)), exp(x), b(x^7, x^8)}:

d := select(t -> t::numeric and t <> 0, degree~(L, x));
F := copy(L);
while indets(F minus {x}) <> {} do
  F := op~(F);
end do:
F minus {x} union d
               {-4, 1, 2, 3, 4, 5, 6, 7, 8}


Finaly, if some exponents are not integers and/or L contains terms of the form something*x^exponent you could use this code:

L := {3, 4, x, x^2+a*x^s, q*x^3, a(x^5), x^(-4), b(a(x^6)), exp(x^(-3/2)), b(r*x^7, x^p)}:

E := select(t -> t::numeric and t <> 0, degree~(L, x)):
F := select(t -> hastype(t, `^` ), L):
while F <> {x} do
  F, R := selectremove(has, F, x);
  F    := map(t -> if t::`*` then select(has, [op(t)], x)[] else op(t) end if, F);
  E    := E union R;
end do:
E;
                /                      -3      \ 
               { -4, 1, 2, 3, 5, 6, 7, --, p, s }
                \                      2       / 

Get_All_Exponents.mw

a := 3:

verify(a, 2..3, 'interval');                    # default means (2, 3)
verify(a, 2..3, 'interval'('closed'));          # here the interval is [2, 3]
verify(a, 2..3, 'interval'('closed'..'open'));  # interval is [2, 3)
verify(a, 2..3, 'interval'('open'..'closed'));  # interval is (2, 3]
                             false
                              true
                             false
                              true

 


You have to fix this big problem before

# Order of eq1 
D_terms  := select(has, indets(eq1, function), D):
Eq1Order := max(degree~(eval(eval(%, f = (eta -> exp(k*eta))), eta=0))); 
                               3

# Total number of IC/BC you impose
ficbc := select(has, {ics, bcs}, f):
NumberIcBc := numelems(%);
                               4

The number of IC/BC is not equal to the order of eq1.
Same thing happen to eq2: a 2nd order ode with 3 IC/BC.

Once you have fix these two problems you will be able (or maybe not...) to solve your differential system:

# Example 1 : eq1 with 3 ICs and eq2 with 2 ICs
# 
# As you then have "free" parameters [a, n, m] use option "parameters".

sol_1 := dsolve(
           {eq1, eq2} union select(has, {ics}, {f, theta})
           , numeric
           , parameters=[a, b, m]
         )
proc(x_rkf45)  ...  end;


# Example 2 : eq1 with 2 ICs + 1 BC and eq2 with 2 ICs.
#
# In this case you solve a coupled bv problem and you cannot use the option "parameters".
# You must instanciate [a, b, m] BEFORE calling dsolve, for instance:

data := [a=1, b=2, m=3]:

sys   := eval({eq1, eq2} union {f(0)=0, D(f)(0)=0, D(f)(10)=1, theta(0) = 1, (D(theta))(0) = b}, data):
sol_2 := dsolve(
           sys
           , numeric
         )

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

This error is quite common when you solve bv problems.
There exist several ways to get rid of this error but they are problem-dependent, in particular they may depend here on the values in data. So i will not waste tume to solve a problem wich could not be the one you are interested in.


Last point: writting things like Cp[f] := 4179 (for instance) is not safe (f is also the name of the function ). Look at this

restart
Cp[f] := 4179
                              4179
f := x^2;
                                2
                               x 
Cp[f]
                                 2 
                             Cp[x ]

Write Cp__f := 4179 instead (same thing for all the quantities whose names nontain [f])

@nm 

r := [[5,7],[],[1,3],[],[5,4]];
remove(type, r, [])


F(z) is not convergent if |z| > 1
 

restart: 
alias(po=pochhammer):
F := z -> sum(po(1,n)^2/po(1/2,n)^2*z^n/(n+1),n=1..infinity):
F(1/2)
                      /              [3  3   ]  1\
             hypergeom|[1, 2, 2, 2], [-, -, 3], -|
                      \              [2  2   ]  2/
F(2)
                            infinity
Digits:=10:
eval(F(z), z=2);
evalf[5](%);
                       /              [3  3   ]   \
            4 hypergeom|[1, 2, 2, 2], [-, -, 3], 2|
                       \              [2  2   ]   /
                       -3.8666 + 3.5627 I

As serie F is not convergent at point z=2 you cannot substutite the value z=2 in the general expression of G(z) because this expression is only valid only if G converges at this point.

A simpler example

n := 'n':
G := z -> sum(z^n, n=0..infinity);
     infinity  
      -----    
       \       
        )     n
z ->   /     z 
      -----    
      n = 0    
'G(2)' = G(2), 'G(z)' = G(z), 'eval(G(z), z=2)' = eval(G(z), z=2);
                                     1                           
   G(2) = infinity,       G(z) = - ------,       eval(G(z), z = 2) = -1
                                   -1 + z                        

 




I found it interesting to show one can find the result without using a while-loop, even though this code is pretty crude and can probably be improved.
 

NT := proc(p, q)
  local ip, iq, ip1, iq1, div, pf, qf, rf:
  ip  := ifactors(p)[2];
  iq  := ifactors(q)[2];
  ip1 := map2(op, 1, ip):
  iq1 := map2(op, 1, iq):
  if `subset`({iq1[]}, {ip1[]}) then 
    div := `intersect`({ip1[]} , {iq1[]});
    pf  := select((x -> op(1, x) in div), ip);
    qf  := select((x -> op(1, x) in div), iq);
    rf  := min(iquo~(map2(op, 2, pf), map2(op, 2, qf)));
    printf("%d divides %d  %d times \n", q, p, rf)
  else 
    printf("%d does not divide %d \n", q, p)
  end if;
end proc:

NT(294912, 8);
8 divides 294912  5 times 

NT(294912, 8*3);
24 divides 294912  2 times 

NT(90,3);
3 divides 90  2 times 

 


 

restart

ineq := log[2](7/10*x) - log[3](3*x-1) <= 0;

ln((7/10)*x)/ln(2)-ln(3*x-1)/ln(3) <= 0

(1)

lineq := lhs(ineq);

plot(
  lineq
  , x=1/3..20
  , color=blue
  , legend=typeset(lineq)
)

ln((7/10)*x)/ln(2)-ln(3*x-1)/ln(3)

 

 

# find values of x such that lineq = 0

s  := solve(lineq);
r  := allvalues(s);
rf := evalf(r);

(1/3)*3^(RootOf(30*exp(_Z)-7*exp(_Z*ln(3)/ln(2))-7)/ln(2))+1/3

 

(1/3)*3^(RootOf(30*exp(_Z)-7*exp(_Z*ln(3)/ln(2))-7, -1.342818283)/ln(2))+1/3, (1/3)*3^(RootOf(30*exp(_Z)-7*exp(_Z*ln(3)/ln(2))-7, 2.453168519)/ln(2))+1/3

 

.3730125034, 16.60731835

(2)

# lineq being continuous for x > 1/3, there is a point x=M between r[1] and r[2]
# where lineq reaches its extremal value.

diff(lhs(ineq), x);
M := solve(%=0);
verify(M, r[1]..r[2], interval);

1/(x*ln(2))-3/((3*x-1)*ln(3))

 

-(1/3)*ln(3)/(ln(2)-ln(3))

 

true

(3)

# At x = M the second ordre derivative of lineq being positive, lineq reaches its minimum:

diff(lineq, x$2);
simplify(eval(%, x=M));
is(%, positive)

-1/(x^2*ln(2))+9/((3*x-1)^2*ln(3))

 

-9*(ln(2)-ln(3))^3/(ln(3)^2*ln(2)^2)

 

true

(4)

# Then lineq < 0 for each x in (r[1], r[2]).
# Or, if you prefer:

x in RealRange(Open(rf[1]), Open(rf[2])) implies lineq < 0;

x in RealRange(rf[1], rf[2]) implies lineq <= 0

`in`(x, RealRange(Open(.3730125034), Open(16.60731835))) implies ln((7/10)*x)/ln(2)-ln(3*x-1)/ln(3) < 0

 

`in`(x, RealRange(.3730125034, 16.60731835)) implies ln((7/10)*x)/ln(2)-ln(3*x-1)/ln(3) <= 0

(5)

 


 

Download By_hand.mw


Begin to look to functions like simplify, normal, factor, collect to cite a few and try to do some things by yourself.
You should find thatthe answer to your first question is obvious.

The one to your second question is not and it is impossible that Maple, or any other CAS, can provide the "simplified" expression you want, just because this simplified form is subjective (meaning it is who claims it is a simple form but others could say another form is simpler).
If guide Maple to the result you want you will get it, but Maple cannot finf it by itself

expr1 := 4*(theta - 1/4)*gamma1^2:
simplify(expr1)
                                         2
                     (4 theta - 1) gamma1 
expr2 := (alpha + gamma1)^2*epsilon^2 + 2*(alpha + gamma1)*(beta + gamma1)*epsilon - 4*alpha*gamma1*theta + beta^2 + 2*beta*gamma1 - 4*(theta - 1/4)*gamma1^2:

Student:-Precalculus:-CompleteSquare(expr2, beta):
op(1, %) + factor(simplify(%-op(1, %)));
                                                    2
          (beta + (alpha + gamma1) epsilon + gamma1)   - 4 gamma1 theta (alpha + gamma1)

solutions.mw

@Kitonum 

parms:=[colour=red,thickness=`if`(leader=0, 3, 0),linestyle =dash]:


Personnaly I prefer writting my frequency table (which is not that difficult) instead of using the build-in procedure which I feel quite limited (IMO).

Here is an example of a procedure (it's name F reflects my lack of imagination):
MyFrequencyTable.mw

It acts differently depending wether you want to consider the data are observations of a discrete process or of a continuous one.
You can also define your own header names (for instance using french words)  manage some colors and chose the display precision for percentages.
It should be easy for you to customize it as you want.

You wrote

p2b==g2b(x-c*t)

instead of

p2b=g2b(x-c*t)

Once corrected you will get what you want.

In case you find that the expressions you get are a little bit obscure or not synthetic enough, this file can give you some ideas for alternative displays forms.mw

As a table is a list of equalities indices=entries, it all depends on what you want these indices and entries are.

For instance it could be natural (@Rouben Rostamian  answer) to think that A is the index of  `A=70` and 70 its entry. But what about "DOT=106"?
Is DOT the corresponding index or is it "DOT"?
Is "106" the corresponding entry or is it 106"?

The attached file presents five different ways (among many others while you do not define indices and entries)  to construct a table from your list st.

restart:

with(StringTools):

st := [`A=70`, `B=17`, `C=27`, `D=37`, `E=74`, `F=57`, `G=67`, `H=08`, `I=81`, `J=28`, `K=38`, `L=48`, `M=58`, `N=68`, `O=90`, `P=19`, `Q=29`, `R=39`, `S=49`, `T=59`, `U=96`, `V=010`, `W=110`, `X=210`, `Y=310`, `Z=410`, "SPACE=105", "DOT=106"]

[`A=70`, `B=17`, `C=27`, `D=37`, `E=74`, `F=57`, `G=67`, `H=08`, `I=81`, `J=28`, `K=38`, `L=48`, `M=58`, `N=68`, `O=90`, `P=19`, `Q=29`, `R=39`, `S=49`, `T=59`, `U=96`, `V=010`, `W=110`, `X=210`, `Y=310`, `Z=410`, "SPACE=105", "DOT=106"]

(1)

sst := convert~(st, string):

 

Option 1

map(x -> Split(x,"="), sst):
map(x -> x[1]=x[2], %):
convert(%, table)

table( [( "M" ) = "58", ( "D" ) = "37", ( "C" ) = "27", ( "V" ) = "010", ( "I" ) = "81", ( "SPACE" ) = "105", ( "R" ) = "39", ( "E" ) = "74", ( "A" ) = "70", ( "L" ) = "48", ( "W" ) = "110", ( "H" ) = "08", ( "S" ) = "49", ( "Z" ) = "410", ( "T" ) = "59", ( "O" ) = "90", ( "F" ) = "57", ( "B" ) = "17", ( "Y" ) = "310", ( "P" ) = "19", ( "K" ) = "38", ( "U" ) = "96", ( "DOT" ) = "106", ( "G" ) = "67", ( "N" ) = "68", ( "Q" ) = "29", ( "X" ) = "210", ( "J" ) = "28" ] )

(2)


Option 2

map(x -> Split(x,"="), sst):
map(x -> x[1]=parse(x[2]), %):
convert(%, table)

table( [( "M" ) = 58, ( "D" ) = 37, ( "C" ) = 27, ( "V" ) = 10, ( "I" ) = 81, ( "SPACE" ) = 105, ( "R" ) = 39, ( "E" ) = 74, ( "A" ) = 70, ( "L" ) = 48, ( "W" ) = 110, ( "H" ) = 8, ( "S" ) = 49, ( "Z" ) = 410, ( "T" ) = 59, ( "O" ) = 90, ( "F" ) = 57, ( "B" ) = 17, ( "Y" ) = 310, ( "P" ) = 19, ( "K" ) = 38, ( "U" ) = 96, ( "DOT" ) = 106, ( "G" ) = 67, ( "N" ) = 68, ( "Q" ) = 29, ( "X" ) = 210, ( "J" ) = 28 ] )

(3)


Option 3

convert(parse~(st), table)

table( [( Q ) = 29, ( Y ) = 310, ( C ) = 27, ( E ) = 74, ( M ) = 58, ( Z ) = 410, ( P ) = 19, ( N ) = 68, ( B ) = 17, ( S ) = 49, ( SPACE ) = 105, ( A ) = 70, ( J ) = 28, ( DOT ) = 106, ( R ) = 39, ( T ) = 59, ( K ) = 38, ( L ) = 48, ( G ) = 67, ( I ) = 81, ( V ) = 10, ( H ) = 8, ( O ) = 90, ( U ) = 96, ( D ) = 37, ( X ) = 210, ( F ) = 57, ( W ) = 110 ] )

(4)


Option 4

map(x -> if x::string then `=`(op(Split(x, "="))) else parse(x) end if, st):
convert(%, table)

table( [( Q ) = 29, ( "SPACE" ) = "105", ( Y ) = 310, ( C ) = 27, ( E ) = 74, ( M ) = 58, ( Z ) = 410, ( P ) = 19, ( N ) = 68, ( B ) = 17, ( S ) = 49, ( "DOT" ) = "106", ( A ) = 70, ( J ) = 28, ( R ) = 39, ( T ) = 59, ( K ) = 38, ( L ) = 48, ( G ) = 67, ( I ) = 81, ( V ) = 10, ( H ) = 8, ( O ) = 90, ( U ) = 96, ( D ) = 37, ( X ) = 210, ( F ) = 57, ( W ) = 110 ] )

(5)


Option 5

map(x -> if x::string then op(1, Split(x, "="))=parse(op(2, Split(x, "="))) else parse(x) end if, st):
convert(%, table)

table( [( Q ) = 29, ( "SPACE" ) = 105, ( Y ) = 310, ( C ) = 27, ( E ) = 74, ( M ) = 58, ( Z ) = 410, ( P ) = 19, ( N ) = 68, ( B ) = 17, ( S ) = 49, ( "DOT" ) = 106, ( A ) = 70, ( J ) = 28, ( R ) = 39, ( T ) = 59, ( K ) = 38, ( L ) = 48, ( G ) = 67, ( I ) = 81, ( V ) = 10, ( H ) = 8, ( O ) = 90, ( U ) = 96, ( D ) = 37, ( X ) = 210, ( F ) = 57, ( W ) = 110 ] )

(6)


....

Download lst.mw

1 2 3 4 5 6 Page 1 of 6