restart: with(Statistics): local E; # Let X some ABSTRACT Random Variable. # # Remark : Yet I'm not able to declare that X is simply a random variable, # that is a random variable without any specified distribution # # So I use X as a DEFINED Random Variable (here Normal) X := RandomVariable(Normal(m,s)); SSNfUkc2JSUqcHJvdGVjdGVkR0kvUmFuZG9tVmFyaWFibGVHNiJJOV9Qcm9iYWJpbGl0eURpc3RyaWJ1dGlvbkdGJg== # Are X and X^2 RandomVariables : type(X, RandomVariable); type(X^2, RandomVariable); SSV0cnVlRyUqcHJvdGVjdGVkRw== SSZmYWxzZUclKnByb3RlY3RlZEc= # I want X^n to be of type RandomVariable for any non negative integer n and, # if U is another RandomVariable, that X*U be of type RandomVariable too. # # So I define a new type 'RV' for "RandomVariables" TypeTools[AddType]( RV, t->`if`(t::`^`, evalb((op(1,t)::Statistics:-RandomVariable) and (op(2,t)::posint)), `if`(t::`*`, evalb((op(1,t)::Statistics:-RandomVariable) and (op(2,t)::Statistics:-RandomVariable)), evalb(t::Statistics:-RandomVariable) ) ) ); type(X, 'RV'); type(X^2, 'RV'); U := RandomVariable(Uniform(p,q)): type(X*U, 'RV'); SSV0cnVlRyUqcHJvdGVjdGVkRw== SSV0cnVlRyUqcHJvdGVjdGVkRw== SSV0cnVlRyUqcHJvdGVjdGVkRw== # Let Y a new random variable defined by Y = f(X); # (f is a deterministic function) Y := X -> f(X); Zio2I0kiWEc2IkYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlLUkiZkdGJTYjOSRGJUYlRiU= # Compute the Taylor expansion TY of Y at "some point" X=a # Note there is an ambiguity in this very common notation because one should better write : # " TY is the Taylor expansion of Y given the value of X is a" TY := convert(taylor(Y(X), X=a, 2), polynom); LCYtSSJmRzYiNiNJImFHRiUiIiIqJi0tSSJERzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YlNiNGJEYmRigsJkkjX1JHNiVGLkkvUmFuZG9tVmFyaWFibGVHRiVJOV9Qcm9iYWJpbGl0eURpc3RyaWJ1dGlvbkdGJUYoRichIiJGKEYo # Define now the operator "MathematicalExpectation" (for short here "E", thus the declaration "local E" above). # # 1) I want E to be linear with respect to objects of type 'RV' : # Example : if Z := a*X+b*Y where X and Y are of type 'RV' and a and b anything, # E(Z) = a*E(X)+b*E(Y) # # 2) When applied to Z = a*X*Y, I want that E(Z) = a*E(X*Y) # # 3) When applied to Z = a*X^n, I want that E(Z) = a*E(X^n) # # I failed in using correctly options 'linear' and 'orderless' E := 'E': define( E, #### 'orderless', E( (a::function)+(b::anything)) = a + E(b), E( (a::`^`)+(b::anything)) = a + E(b), E( (a::symbol)+(b::'RV')) = a + E(b), E( (b::'RV')-(a::symbol)) = E(b)-a, E( (a::function)*(b::anything)) = a * E(b), E( (a::`^`)*(b::anything)) = a * E(b), E( (a::symbol)*(b::'RV')) = a * E(b) ); # Let me apply E to TY E(TY); LCYtSSJmRzYiNiNJImFHRiUiIiIqJi0tSSJERzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YlNiNGJEYmRigsJi1JIkVHNiRJJ19sb2NhbEdGLi9GNUktX240Mzc0NjI3MjY0R0YlNiNJI19SRzYlRi5JL1JhbmRvbVZhcmlhYmxlR0YlSTlfUHJvYmFiaWxpdHlEaXN0cmlidXRpb25HRiVGKEYnISIiRihGKA== # Sounds good (it is likely that my definitions of 'RV' and 'E' are neither generic nor robust ...) # # Now apply E to the square of TY # # Its is expected that E(TY^2) will return an expression of the form expr := A + B*E(X) + C*E(X^2) # If one substract to expr the square of E(TY) (which is of the same form) then, for a suitable choice # of a (hre the Mean, or mathematical expectation, of X), then E(TY^2)-(E(TY)^2) is a first order # estimation of the variance of Y # # # More details may be found in : # https://en.wikipedia.org/wiki/Propagation_of_uncertainty # https://www.physics.ohio-state.edu/~gan/teaching/spring04/Chapter4.pdf # # but it is a widely used method for quickly assess values of the mean and variance of a function # of many random variable # TY2 := expand(TY*TY); LC4qJC1JImZHNiI2I0kiYUdGJiIiIyIiIiooRiRGKi0tSSJERzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YmNiNGJUYnRipJI19SRzYlRjBJL1JhbmRvbVZhcmlhYmxlR0YmSTlfUHJvYmFiaWxpdHlEaXN0cmlidXRpb25HRiZGKkYpKihGJEYqRixGKkYoRiohIiMqJkYsRilGM0YpRioqKEYsRilGM0YqRihGKkY4KiZGLEYpRihGKUYq E(TY2); LCYqJC1JImZHNiI2I0kiYUdGJiIiIyIiIi1JIkVHNiRJJ19sb2NhbEclKnByb3RlY3RlZEcvRi5JLV9uNDM3NDYyNzI2NEdGJjYjLCwqJi0tSSJERzYkRi9JKF9zeXNsaWJHRiY2I0YlRidGKUkjX1JHNiVGL0kvUmFuZG9tVmFyaWFibGVHRiZJOV9Qcm9iYWJpbGl0eURpc3RyaWJ1dGlvbkdGJkYpRioqJkY1RilGKEYpRioqKEYkRipGNUYqRjtGKkYpKihGJEYqRjVGKkYoRiohIiMqKEY1RilGO0YqRihGKkZCRio= # ... you see here that E does not seem to operate. # probably its construction is too simple to account for all the situations, or there is some precedence # order between the type evaluation ('RV') and the operator 'E' itself ?