8 years, 69 days

## @fabs  Here is the edited file whi...

@fabs

Here is the edited file which includes the case X=x1*tauderivatives.mw
In fact the trick works whatever the function of x1 and tau.

## @Carl Love  Thanks Carl, I underst...

@Carl Love

Thanks Carl, I understand better now.

## @Carl Love  Than you Carl. This is...

@Carl Love

Than you Carl.
This is almost what I did just after havng read @vv's answer.

```restart
phi := x -> (x^(-theta)-1)/theta:
ihp := u -> solve(phi(x)=u, x):
invfunc[phi] := ihp:

C := (u, v) -> (phi@@(-1))(phi(u)+phi(v)):
simplify(C(u, v)) assuming theta >= -1, theta <> 0;
/    1  \
|- -----|
\  theta/
/      (-theta)    (-theta)\
\-1 + u         + v        /
```

The only slight difference (but maybe of great importance???) is that I did not use unprotect(invfunc).
I only checked out that this table indeed contained phi=ihp.

## @vv  Thank you vv. I'm feeling...

@vv

Thank you vv.
I'm feeling sorry for being so stupid because I had looked at the invfunc help page but obviously not carefully enough.

## @acer Thank you for this suggestion...

Thank you for this suggestion and the related comments

I only knew about the operator-form calling sequence for plotting commands ot in the case of Optimization, where I use it on a regular basis.
I didn't know it could also be used int too.

And yes "the approach can break due to some overlooked extra evaluation which strips off a pair of such quotes": I unfortunately made this bad experience a lot of times.

## @Preben Alsholm  Thanks you Preben...

@Preben Alsholm

Thanks you Preben.

## Thanks...

@Preben Alsholm

Sometimes I write my procedures the way you suggest, in particular when these procedures are used in a plot command.
But I also noted that these definitions both gave the same correct result:

```MyProc1 := proc(u)
if not u::realcons then return 'procname'(_passed) end if;
fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = u, x)
end proc:

plot(MyProc1(u), u=0.1..0.9):

MyProc2 := proc(u)
fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = u, x)
end proc:

plot('MyProc2'(u), u=0.1..0.9):

```

Son lazzy as I am, I often use the second definition.
This is why I used it in my question.

So, it seems that defining the procedure as in MyProc1 is "safer' than in MyProc2. Am I right?
Do you systematically recommend proceeding as in MyProc1 ?

By the way: the computational time comes from the fact that

`fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = q, x)`

is equal to -infinity if q=0 and +infinity if q=1. In practice I would integreate within the range 0.01..0.99 for instance.

For the record here is the method I used to use to numerically evaluate J2 before I thought of using evalf/Int instead (last block in the attached file).
I suppose that a judicious choice of the method in evalf/Int can provide the result fastly then the default method used in J3 ?

 > restart:
 > J2 := proc()   local z:    z := proc(q1, q2)        local x1, x2;      if not [q1,q2]::list(realcons) then return 'procname'(_passed) end if;      x1 := fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = q1, x):      x2 := fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = q2, x):      0.2652582384*exp( 2*x1*x2 - x1^2 - x2^2 )   end proc:   evalf(Int(z(q1, q2), q1=0.01..0.99, q2=0.01..0.99)) end proc: CodeTools:-Usage( J2() );
 > J3 := proc()   local z:    z := proc(q1, q2)        local x1, x2;      uses Statistics:      if not [q1,q2]::list(realcons) then return 'procname'(_passed) end if;      x1 := Quantile(Normal(0, 1), q1, numeric):      x2 := Quantile(Normal(0, 1), q2, numeric):      exp( 2*x1*x2 - x1^2 - x2^2 )   end proc:   evalf(Int(z(q1, q2), q1=0.01..0.99, q2=0.01..0.99)) end proc: CodeTools:-Usage( J3() );
 memory used=0.76GiB, alloc change=44.01MiB, cpu time=14.31s, real time=14.55s, gc time=689.80ms
 (1)

A more efficient way

 > StringTools:-FormatTime("%H:%M:%S"); h := 0.001: P := [seq(0.01..0.99, h)]: Q := CodeTools:-Usage( Statistics:-Quantile~(Normal(0, 1), P, numeric) ): S:= 0: for q1 in Q do   for q2 in Q do     S := S + exp( 2*q1*q2 - q1^2 - q2^2 )   end do: end do: S * h^2; StringTools:-FormatTime("%H:%M:%S");
 memory used=16.68MiB, alloc change=0 bytes, cpu time=291.00ms, real time=292.00ms, gc time=0ns
 (2)
 >

## Optional argument...

@Ronan

Writing

`{clr::string:="Blue"}`

in the argument list of a procedure means that clr is to be considered as an optional argument whose type is string and default value "Blue".
If you do not mention clr when you execute the procedure spread2, clr will take inside this procedure the default value "Blue".
If you want the procedure uses "Red" instead of "Blue", you have to set ecplivitely clr="Red" when you call it.

Optional arguments are intensively used in Maple: you may see them as an annoyance for they impose writing clr="Red" instead of simply "Red", but their advantage is that the order ofthe optional parameters doesn't matter (so you can write indifferently clr="Red", prnt=false, or prnt=false, clr="Red").
I greatly appreciate this feature, but it remains a personal opinion.

## @acer  Thank you acer. Your eval...

@acer

Thank you acer.
Your

`eval(y, g=''f'')`

suits me perfectly.

By the way I sent at @Preben Alsholm a file where the real problem is described, in case you would be interested.

## @Preben Alsholm  Thanks four answe...

@Preben Alsholm

"why not simply start with f=F?" simply because this example was notional (maybe oversimplified).
Here is a more detailed example which is closer to the reality (note that writing it I realized there was a simple way [end of the file] to answer my own question!)

 > restart
 > with(LinearAlgebra): with(Statistics):
 > phi := unapply(CDF(Normal(0, 1), x), x);
 (1)
 > Sigma := Matrix(2\$2, [1, rho, rho, 1]): X := Vector(2, symbol=x):
 > phi__Sigma := exp(-1/2*(X^+ . Sigma^(-1) . X)) / (2*Pi*sqrt(Determinant(Sigma)));
 (2)
 > # Let psi the CDF of some continuous random variable V. # We define C y C__psi := eval(phi__Sigma, [seq(x[i]=((phi@@(-1))@cdf)(q[i]), i=1..2)]);
 (3)
 > # In the specific case V is a standard gaussian random variable we have cdf(V) = phi C__psi := eval(phi__Sigma, [seq(x[i]=((phi@@(-1))@phi)(q[i]), i=1..2)]);
 (4)
 > # But cdf := x -> phi(x): C__psi := eval(phi__Sigma, [seq(x[i]=((phi@@(-1))@cdf)(q[i]), i=1..2)]);
 (5)

Writing this more realistic example revealed me how to proceed in a simpler way

 > C:= psi ->  eval(phi__Sigma, [seq(x[i]=((phi@@(-1))@psi)(q[i]), i=1..2)]):
 > C(psi); C(phi);
 (6)
 >

## dl1dsd, dl1dsv, dl1dg and a few more t...

 > restart
 > local Gamma, gamma:
 > Gamma := phi(sigma__d): F := f(Gamma)*sigma__v/sigma__d: dl1dsd := convert(diff(F, sigma__d), diff): dl1dsd := subs(diff(phi(sigma__d), sigma__d)=freeze(diff(phi(sigma__d), sigma__d)), dl1dsd): Gamma := 'Gamma': dl1dsd := eval(dl1dsd, phi = (sigma__d -> Gamma)): dl1dsd := thaw(dl1dsd): dl1dsd := eval(dl1dsd, phi = (sigma__d -> gamma*sigma__v*sigma__d));
 (1)
 > Gamma := phi(sigma__v): F := f(Gamma)*sigma__v/sigma__d: dl1dsv := convert(diff(F, sigma__v), diff): dl1dsv := subs(diff(phi(sigma__v), sigma__v)=freeze(diff(phi(sigma__v), sigma__v)), dl1dsv): Gamma := 'Gamma': dl1dsv := eval(dl1dsv, phi = (sigma__v -> Gamma)): dl1dsv := thaw(dl1dsv):   dl1dsv := eval(dl1dsv, phi = (sigma__v -> gamma*sigma__v*sigma__d));
 (2)
 > local gamma: Gamma := phi(gamma): F := f(Gamma)*sigma__v/sigma__d: dl1dg := convert(diff(F, gamma), diff): dl1dg := subs(diff(phi(gamma), gamma)=freeze(diff(phi(gamma), gamma)), dl1dg): Gamma := 'Gamma': dl1dg := eval(dl1dg, phi = (gamma -> Gamma)): dl1dg := thaw(dl1dg): dl1dg := eval(dl1dg, phi = (gamma -> gamma*sigma__v*sigma__d));
 (3)
 >

If you want to get more synthetic displays for der1 and der2  here is a hint:

 > Lambda__1 := RootOf(8*_Z^4 + 12*Gamma*_Z^3 + (5*Gamma^2 - 4)*_Z^2 - 4*Gamma*_Z - Gamma^2):
 > # A more synthetic representation of der1: der1 := diff(Lambda__1,Gamma): Diff('Lambda__1',Gamma) = collect~(normal(eval(der1, Lambda__1 = Lambda[1])), Lambda[1]);
 (1)
 > # A more synthetic representation of der2: der2 := diff(der1,Gamma): Diff('Lambda__1', Gamma\$2) = collect~(normal(eval(der2, Lambda__1 = Lambda[1])), Lambda[1]);
 (2)
 >

Last point: lambda__1 has 4 roots (4 real, or 2 real and 2 omplex conjugate, or 2 complex conjugate and 2 other complex conjugate too).
When you do plot(lambda__1, Gama=1..10) you plot only one of them:

```eval(Lambda__1, Gamma=0):
evalf([allvalues(%)]);

[0., 0., 0.7071067810, -0.7071067810]

eval(Lambda__1, Gamma=3):
evalf([allvalues(%)]);

[0.4958019689, -2.302318663 + 0.7071675271 I, -0.3911646429,
-2.302318663 - 0.7071675271 I]

eval(Lambda__1, Gamma=10):
evalf([allvalues(%)])

[0.4628298392, -7.515971909 + 2.487925010 I, -0.4308860212,
-7.515971909 - 2.487925010 I]
```

## Why Quantity?...

@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;
 > 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;
 (1)
 >

## Statistics and Units...

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')
 (1)
 > Mean(Y);
 (2)
 > Variance(Y);
 (3)
 > Kurtosis(Y);  # should be dimensionless
 (4)
 > restart
 > with(Statistics):
 > Y := RandomVariable(Normal(mu, sigma))*Units:-Unit('m')
 (5)
 > Mean(Y);
 (6)
 > Variance(Y);  # units bad placed
 (7)
 > Kurtosis(Y);
 (8)
 >