mmcdara

7891 Reputation

22 Badges

9 years, 54 days

MaplePrimes Activity


These are replies submitted by mmcdara

@Nicole Sharp 

If you want to use ScientificErrorAnalysys to get the mean value and the error on rho[cubewano] you must do the things as you did.
See this notional example:
 

restart:

use ScientificErrorAnalysis, Units in
  X := Quantity(10, 2)*Unit(m);
  Y := Quantity( 3, 4)*Unit(kg);
  Z := Quantity(10, 1)*Unit(s);

  F := Y*X/Z^2:

  E    := combine(F, 'errors'):
  mean := GetValue(E);
  std  := GetError(E);
end use;

ScientificErrorAnalysis:-Quantity(10, 2)*Units:-Unit('m')

 

ScientificErrorAnalysis:-Quantity(3, 4)*Units:-Unit('kg')

 

ScientificErrorAnalysis:-Quantity(10, 1)*Units:-Unit('s')

 

ScientificErrorAnalysis:-Quantity(3, 4)*Units:-Unit('kg')*ScientificErrorAnalysis:-Quantity(10, 2)*Units:-Unit('m')/(ScientificErrorAnalysis:-Quantity(10, 1)^2*Units:-Unit('s')^2)

 

ScientificErrorAnalysis:-Quantity((3/10)/Units:-Unit('s')^2, (1/50)*418^(1/2)/Units:-Unit('s')^2)*Units:-Unit('kg')*Units:-Unit('m')

 

Error, (in ScientificErrorAnalysis:-GetValue) expecting a quantity-with-error structure, but got ScientificErrorAnalysis:-Quantity((3/10)/Units:-Unit(s)^2, (1/50)*418^(1/2)/Units:-Unit(s)^2)*Units:-Unit(kg)*Units:-Unit(m)

 

use ScientificErrorAnalysis, Units in
  X := Quantity(10*Unit(m) , 2*Unit(m) );
  Y := Quantity( 3*Unit(kg), 4*Unit(kg));
  Z := Quantity(10*Unit(s) , 1*Unit(s) );

  F := Y*X/Z^2:

  E    := combine(F, 'errors'):
  mean := GetValue(E);
  std  := GetError(E);
end use;

ScientificErrorAnalysis:-Quantity(10*Units:-Unit('m'), 2*Units:-Unit('m'))

 

ScientificErrorAnalysis:-Quantity(3*Units:-Unit('kg'), 4*Units:-Unit('kg'))

 

ScientificErrorAnalysis:-Quantity(10*Units:-Unit('s'), Units:-Unit('s'))

 

ScientificErrorAnalysis:-Quantity(3*Units:-Unit('kg'), 4*Units:-Unit('kg'))*ScientificErrorAnalysis:-Quantity(10*Units:-Unit('m'), 2*Units:-Unit('m'))/ScientificErrorAnalysis:-Quantity(10*Units:-Unit('s'), Units:-Unit('s'))^2

 

ScientificErrorAnalysis:-Quantity((3/10)*Units:-Unit('kg')*Units:-Unit('m')/Units:-Unit('s')^2, (1/50)*418^(1/2)*Units:-Unit('kg')*Units:-Unit('m')/Units:-Unit('s')^2)

 

(3/10)*Units:-Unit('kg')*Units:-Unit('m')/Units:-Unit('s')^2

 

(1/50)*418^(1/2)*Units:-Unit('kg')*Units:-Unit('m')/Units:-Unit('s')^2

(1)

 


 

Download Quantity_with_Unit.mw

@Nicole Sharp 

You < also noticed that "with(Statistics) : Mean(%)" doesn't seem to work with Units either >.

(version used: Maple 2015, but I guess this works also in newer versions)

restart

with(Statistics):

Y := RandomVariable(Normal(2, 3))*Units:-Unit('m')

_R*Units:-Unit('m')

(1)

Mean(Y);

2*Units:-Unit('m')

(2)

Variance(Y);

9*Units:-Unit('m')^2

(3)

Kurtosis(Y);  # should be dimensionless

3

(4)

restart

with(Statistics):

Y := RandomVariable(Normal(mu, sigma))*Units:-Unit('m')

_R*Units:-Unit('m')

(5)

Mean(Y);

mu*Units:-Unit('m')

(6)

Variance(Y);  # units bad placed

Units:-Unit('m')^2*sigma^2

(7)

Kurtosis(Y);

3

(8)

 

Download Statistics_Units.mw

@vv 

I agree. 
I titled my answer "The limit is not a number" because I didn't want to waste time developing this point and explaining why he couldn't get the expected result by using "limit".
Your last point is very pertinent.

Thanks for your comment.


HERE

What exactly do you want more?

@acer 

I understand your "When none of your piecewise conditions holds it will evaluate to its default value, which by default is zero.", but the plot you get by adding undefined in the piecewise definition is wrong.
The function g being exactly equal to 0 on the "big red block", it seems natural, default value or not, that implicitplot (try contourplot instead with option countours=[0]) draws a line to "fill" this rectangle.
Simply give a look to the fie I attached to my answer. 

It is therefore purely fortuitous that the undefined option produces the layout that the OP desires: the graph is wrong and it's exactly this wrong  graph he wanted for aesthetic reasons.
Here is an example to (try to) convince you this option gives wrong results: gnull_2.mw

@Beavis 
I have edited my answer right now

@Kitonum 

Excellent, thank you

@acer 

... and if I had paid more attention I would have observed that 

Mspec := Specialize(B, [p=1/2])*Specialize(X, [mu=-3, sigma=1]) + (1-Specialize(B, [p=1/2]))*Specialize(Y, [nu=3, tau=1]);
                   Mspec := _R2*_R3+(1-_R4)*_R5

clearly indicates that one Specialize(B, [p=1/2]) is _R2 while the other is _R4.

Whatever your "demonstration" is perfectly clear.

If it's not too much to ask I would have two questions:

  1. I use the Statistics package for years and I always wonder  why RandomVariable creates a protected name of rhe form _Rxxx?
    It seems to me that this is quite a unique case in Maple isn't it?
    These _Rxxx names look more like aliases of the names to which RandomVariables are assigned (am I right?) and they make difficult to understand expressions which involve random variables (relation (7) in the attached file)  
     
  2. Second question: why is the type RandomVariable is not inherited?

Thank you for your attention.

restart:

with(Statistics):

X := RandomVariable(Uniform(0, 1))

_R

(1)

equivalence := anames(RandomVariable) = eval(anames(RandomVariable))

X = _R

(2)

Mean(X), Mean(_R)

1/2, 1/2

(3)

PDF(X, x), PDF(_R, x)

piecewise(x < 0, 0, x < 1, 1, 0), piecewise(x < 0, 0, x < 1, 1, 0)

(4)

attributes(X);
attributes(_R);

protected, RandomVariable, _ProbabilityDistribution

 

protected, RandomVariable, _ProbabilityDistribution

(5)

Y := RandomVariable(Uniform(0, 1));

_R0

(6)

Z := X+Y

_R+_R0

(7)

A := _R + _R0

_R+_R0

(8)

Mean~([Z, A])

[1, 1]

(9)

type(Z, RandomVariable);
type(A, RandomVariable);

false

 

false

(10)
 

 

Download RV.mw

@sursumCorda 

Thanks for this detailed information

@vv 

Thanks

@Preben Alsholm 

Thank you Preben.
The algorithm I was looking for corresponds indeed to procedure  `convert/real_rat`.

@dharr 

Thanks dharr.
I actually thought of it myself shortly after posting my question, but didn't take the trouble to correct it.
Nevertheless I vote up.
 

@C_R 

That is indeed a good idea I hadn't thought of (I vote up).

Regarding your comment about units, this is a point I haven't addressed, but it could be of some interest.
I need to think about it more deeply, but I'll keep it in mind. Anyway, thanks for the tip

@Carl Love 

After correction I get these results: horner_stat.mw

Two observations:

  1. The evaluation of the Horner factored polynomial is about twice faster than than the evaluation of the polynomial itself and uses about half the latter's.
    restart
    with(Statistics):
    U := RandomVariable(Uniform(0, 1)):
    P := randpoly(U, dense, degree=500): H := convert(P, horner, U):
    p := unapply(P, U):
    h := unapply(H, U):
    
    SU := Sample(U, 10^4):
    SP := CodeTools:-Usage( p~(SU) ):
    SH := CodeTools:-Usage( h~(SU) ):
    memory used=226.30MiB, alloc change=76.01MiB, cpu time=2.35s, real time=2.35s, gc time=219.37ms
    memory used=152.99MiB, alloc change=32.00MiB, cpu time=1.34s, real time=1.31s, gc time=50.12ms
    So the difference is not that large than the one you display in you answer to @Aliocha
    My question is: Why?
  2. Can you explain me a direct sampling of the two polynomials:
    1. give almost the same performances?
      randomize(170478831571186):
      
      Sh := CodeTools:-Usage( Sample(H, 10^4) ):
      
      randomize(170478831571186):
      
      S  := CodeTools:-Usage( Sample(P, 10^4) ):
      memory used=0.69MiB, alloc change=0 bytes, cpu time=619.00ms, real time=619.00ms, gc time=0ns
      memory used=505.40KiB, alloc change=0 bytes, cpu time=681.00ms, real time=681.00ms, gc time=0ns
      
    2. are much faster than the evaluations of these polynomials on a  U sample (point 1 bove)?

Finally, as the data type emphasize or lessen the differences between Horner and non Horner:

# With Integers

SU := LinearAlgebra:-RandomVector(10^4, generator=-10^3..10^3):
SP := CodeTools:-Usage( p~(SU) ):
SH := CodeTools:-Usage( h~(SU) ):
memory used=4.53GiB, alloc change=64.00MiB, cpu time=8.58s, real time=6.30s, gc time=3.76s
memory used=2.79GiB, alloc change=0 bytes, cpu time=4.30s, real time=2.99s, gc time=2.10s

# With Rationals

S1 := LinearAlgebra:-RandomVector(10^3, generator=1..10^3):
S2 := LinearAlgebra:-RandomVector(10^3, generator=1..10^3):
S3 := 2 *~ LinearAlgebra:-RandomVector(10^3, generator=0..1) -~1:
SU := S3 *~ S1 /~ S2:

SP := CodeTools:-Usage( p~(SU) ):
SH := CodeTools:-Usage( h~(SU) ):
memory used=1.78GiB, alloc change=0 bytes, cpu time=13.81s, real time=13.24s, gc time=925.67ms
memory used=0.95GiB, alloc change=0 bytes, cpu time=1.38s, real time=1.09s, gc time=458.98ms
-

when is Horner factorization useful and when it's useless?

# With Irrationals

S1 := LinearAlgebra:-RandomVector(10^4, generator=1..10^3):
SU := sqrt~(S1):

SP := CodeTools:-Usage( p~(SU) ):
SH := CodeTools:-Usage( h~(SU) ):
memory used=3.88GiB, alloc change=32.00MiB, cpu time=17.81s, real time=11.56s, gc time=8.61s
memory used=1.76GiB, alloc change=32.00MiB, cpu time=12.25s, real time=8.82s, gc time=5.09s

SP[1]; 

3316745734089766397592630110471859650386443147953483621189845470\

  90097763832529787408149010984184954904189345464801382358281506\

  33776102815748690885532435858761523283924372415325847138335927\

  21476776724692220329022997486496244374362690018578817495888580\

  89276047423467309118046638389871685207180192851631794711037526\

  68226650565956387221753112011216542607616820652974153405307803\

  11762352016961376448324782073540088181636152149702927597228541\

  10972110620460287943566861248221758332972512841962698051882950\

  87585311899955236753512717342312862424600797702334851148866019\

  87487905716563591872683982570084096256655551807913554703328334\

                                  (1/2)                          
  2684004425077425823763313576 401      - 4360993989574992282441\

  74655444335506239538092992960307581344624831828308577448248585\

  86643733234979191263532398267659542635477007073264697573690352\

  41184808202886271956377113251918126274819079448214120087602580\

  40976153761457648322734186550766933126965981655745312867389989\

  47273948044111754646777788245934072208799317202380511846729042\

  84664349281405519565375009208274681063255302780948520203565761\

  02283632683447496819722167297580873880176488752482788795280329\

  32859059120794697071801590589198678945797208588083859312946103\

  31831105107339972858373191659906793831837170122360471704338200\

  49578949360502838215949298234767556200734638282783444666508862\

  7426855551

SH[1];
Error, while processing result

@Carl Love 

I've just updated my last reply

First 31 32 33 34 35 36 37 Last Page 33 of 154