epostma

1339 Reputation

16 Badges

11 years, 179 days
Maplesoft

I am a Senior Software Developer in the Math Group, working mostly on the Maple library. I have been working at Maplesoft since 2007, mostly on the Statistics and Units packages and on contract work with industry users. My background is in abstract algebra, in which I completed a PhD at Eindhoven University of Technology. During my studies I was always searching for interesting questions at the crossroads of math and computer science.

When I came to Canada in 2007, I had nothing but a work permit, some contacts at Maplesoft (whom I had met at a computer algebra conference a year earlier), and a plan to travel around beautiful Canada for a few months. Since Maplesoft is a company that solves more interesting math and computer science related questions than just about anywhere else, I was quite eager to join them, and after about three months, I could start.

MaplePrimes Activity


These are answers submitted by epostma

Hi mmcdara,

I'm afraid I can't reproduce your finding that the Student[Statistics] package does this better than the top level Statistics package; in Maple 2018, if I do:

with(Student[Statistics]):
R := NegativeBinomialRandomVariable(1.965, 0.003);
Quantile(R, 0.98);

I get no answer, just like with the top level Statistics package. In fact, that is to be expected: if you run the command showstat(Student[Statistics][Quantile]), you can see that it directly calls the Quantile command in the top level Statistics package. (That's how the Student[Statistics] package is set up in general - it's meant to be a slightly simpler user interface to the same functionality. Among the few bits of functionality that are exclusive to Student[Statistics] are plots associated with quantities (activated with the output option) and the TestsGuide command).

Maybe more importantly, as I mention in my response to Teep's post, I have a fix to this bug in the works.

Erik Postma
Manager, mathematical software group.

Hi Teep,

Thank you for this report; I saw it today and I have a fix in progress. It looks like the biggest problem is that Maple uses code here that is meant for continuous distributions instead of discrete distributions, for doing root finding, trying to find the value where the CDF crosses the value 0.98. The code fairly quickly finds that at 1919.998, the CDF is 0.9799689, and at 1920.42, the CDF is 0.9800204, but then gets stuck trying to find an epsilon-sized interval in between that contains it. It then tries an alternate method that is a very bad idea.

With this fix, Maple will use a discrete interpolation method for discrete distributions. We've had the code for this since the Statistics package was first written, but it wasn't used in this branch of code. The fix should make it into Maple 2019, unless something unexpected happens.

 

Erik Postma
Manager, mathematical software group.

Thanks for this report, _Maxim_. This issue appears to be fixed in our internal development version. I can't promise any particular future version the fix will show up in, but it should make it out eventually!


Thanks for this report, _Maxim_ - it looks like this issue is fixed in Maple 2018.


As a (partial?) workaround, there is

dsolve(diff(f(z), z) - f(z)/z = 0, f(z), series);

which returns  _C1*z*(1 + O(z^6)) with the default setting for Order.

Unfortunately, this is a known limitation in the current versions of Maple. There are various stages of evaluation of a module definition, and at the time the option option object(ParentObject) is being processed, the ParentObject cannot yet be recognized as a valid object. I am not sure we can fix this in the near future, and all the potential workarounds that I'm aware of have significant downsides.

One such workaround is of course to make the objects top-level entities, instead of living inside a module. In this case, you lose the name space separation that a module provides.

Another would be to define the objects after defining the package:

module Package()
  option package;
  export ParentObject, ChildObject;
end module:

unprotect('Package:-ParentObject');
Package:-ParentObject := module()
  option object;
end module:
protect('Package:-ParentObject');

unprotect('Package:-ChildObject');
Package:-ChildObject := module()
  option object(Package:-ParentObject);
end module:
protect('Package:-ChildObject');

This works well if there is nothing else inside the package, but:

  1. It is clumsy with the protect / unprotect pairs.
  2. If there are other members of Package, then you don't have lexical access to them. In particular, locals of Package are invisible (unless you use kernelopts(opaquemodules=false), which has its own problems).
  3. As a consequence, it requires that all classes that you define are exports of Package. Otherwise, you might for example want to have a common parent class of a number of derived classes be a local that's invisible to the outside world.

Erik Postma
Manager, mathematical software group.


Hi tbfl1919,

 

Let's first investigate what exactly is going on.

The functions V__tot and similar all take a unitless quantity t as their argument, as you can see by evaluating, for example, V__tot(1.3). In the integral, however, you specify the lower and upper bounds for t, the argument of V__tot, as T__1 and T__2, which are 0 and 1 ms. (I should note here that 0 is a special case - it can only be unitless in Maple for technical reasons - but the int command interprets the range as expressed in units of time.) So the int command tries to substitute time values for t, which causes the sin command inside V__tot to fail, as it does on the command line if you evaluate V__tot(1.3*Unit(ms)). Maple decides that that's an expression it can't integrate and merely takes the unit, V2, out of the integral. So we have now evaluated the integral to some unevaluated integral times the unit V2. This is then divided by T__2 - T__1, a duration, to obtain an unevaluated integral with the unit V2/s. The square root of this then has the unit V/s1/2, which is equivalent to the unit you see. It could be argued that this is correct, because if we could specify the integrand better, the integral would come out to be a duration (it is a unitless quantity, integrated over a period of time, so the dimension is time), its square root therefore has the unit s1/2, and multiplying that with the unit we see gives the unit m2 kg / s3 / A, which is the same as V.

 

However, this is all of course incredibly clunky and it doesn't lead to Maple giving you the answer you're looking for. So let's talk about solutions. The first option is to make the integral strip off the time unit and add it back on manually. I think the appropriate unit is ms, correct? So we can do that by specifying the lower and upper bounds as T__1/Unit(ms) and T__2/Unit(ms) and multiplying the integral by Unit(ms). This appears to lead to a correct result (of sqrt(2)/2 V). But this is, in some sense, cheating: we just remove the units manually and then put them back on. It's possible to do better.

 

The second option is to make V__tot accept an expression in terms of time, rather than a unitless quantity. For this we need to do two things:

  1. Redefine V__tot_AC(t) := A * sin(omega * t);
  2. Be very careful to never submit a unitless quantity into V__tot or V__tot_AC.

The latter issue appears, for example, when we try to plot these functions: the only way I can think to do this is to change the command to plot(V__tot(t * Unit(ms)), t = 0 .. 5) (and similar for V__tot_AC). It also appears when we want to do the integration: just entering V__tot(t) already causes an error. However, we can prevent the error by using unevaluation quotes (single forward quotes) around V__tot(t): int('V__tot(t)'2, t = T__1 .. T__2) evaluates no problem. This leads to the same result of sqrt(2)/2 V.

 

You don't write what version of Maple you are using; if it is version 2017, you can use another trick to get around some of the problems here. In 2017 we introduced the Units[Simple] package as an alternative to Units[Standard]. The most important difference between the two is the following: the Units[Standard] package assumes that every unassigned variable, such as t, represents a unitless quantity. The Units[Simple] package does not. This restriction may have played a role in your decision to define V__tot_AC(t) the way you did, where you included the unit for t in the definition, rather than it just being part of the value t that the user submits. Indeed, with the Units[Standard] package, you cannot evaluate sin(omega * t), if t is unassigned and omega has a unit of frequency, because the argument to the sine function is then assumed to have a unit of frequency, which is illegal. On the other hand, if you load the Units[Simple] package, it will see the unassigned variable t and not assume anything about it. (It does assume that every unassigned variable has a single well-defined dimension within one expression; if you evaluate (t + 5 kg) * (t + 5 ms), then Maple will complain: the first factor implies that t must be a mass and the second that t must be a duration.) So compared with the Units[Standard] solution, you don't need to worry about what arguments you submit to V__tot so much: if the argument can be a duration for some assignment of units to your variables, then Maple will accept it. So item 2. above is not so important anymore. In particular, you get the correct answer for V__tot_rms by only removing the unit ms from the definition of V__tot_AC, no quoting needed anywhere.

 

Finally, for the upcoming version of Maple (which is still some time away!) we have overhauled plotting with units somewhat; if you use my suggested definition of V__tot and you use Units[Simple] so that you don't need to worry about evaluating V__tot on symbolic arguments, you will be able to get the correct plot just by writing plot(V__tot(t), t = 0*Unit(ms) .. 5*Unit(ms)) or, if you prefer, plot(V__tot, 0*Unit(ms) .. 5*Unit(ms)). We feel this is a more natural way to use units in plots.

If you replace the last loop with:

e2X := 0;
for a from 1 to 1 do
  for alpha from 0 to k-1 do
    for b from 0 to k-1 do
      for beta from 0 to k-1 do
        f := 2*S*(d[beta+1,alpha+1]*W(beta))*L(alpha)*d[a+1,b+1]*L(b);
        e2X := e2X - f;
      end do;
    end do;
  end do;
end do:
e2 := int(e2X, z = -h/2..h/2);

then you save about half of the running time.

Kitonum has a good suggestion in their answer. Here is an alternative: the typematch command. You pass it the expression to be tested, and a pattern. This pattern consists of type names with a double colon and a name. For example, to match a product of two sums, each with two terms, you could match to the expression (a::algebraic + b::algebraic) * (c::algebraic + d::algebraic). The command returns true if and only if the expression matches the pattern, and in that case it assigns the matching values to the variables occurring in the pattern; in this case to a,b,c,d. See the documentation here: https://www.maplesoft.com/support/help/Maple/view.aspx?path=typematch. One thing that's explained there is how to deal with evaluation problems if the variables in your pattern are assigned values, e.g. if you've used the same pattern successfully before. (It involves unevaluation quotes.)

 

If you don't care for the matching values, just for the yes-or-no result, you can use the type command. It takes a structured type as its second argument, which is basically like the second argument of typematch without the variables. A tricky thing here is that Maple automatically simplifies something like algebraic+algebraic to 2*algebraic, and you don't want that to happen with your pattern.

 

In either case, you may need to study structured types a bit too come up with the right pattern/type. They are explained (rather succinctly) on these help pages: http://www.maplesoft.com/support/help/Maple/view.aspx?path=type/structure and http://www.maplesoft.com/support/help/Maple/view.aspx?path=type. Let us know if you have trouble coming up with the right type to check!

What's nice about the probability distribution is that it can be factored into three factors, each dependent only on x, y, or z, respectively: f = exp(-x^2) * exp(-y^2) * exp(-z^4). That means you can sample the x, y, and z-values independently.

For x and y, it's quite easy: they are normal distributions (centered at 0). We just need to determine the standard deviation - which would be sqrt(2)/2, as you can verify as follows:

with(Statistics):
PDF(Normal(0, sqrt(2)/2), x);
# returns: exp(-x^2)/sqrt(Pi)

For z, we need to expend a little more effort. You need to define a custom distribution, and for this you need to find out the scaling constant needed to make the density integrate to 1. We compute:

c := int(exp(-z^4), z=-infinity .. infinity);
# returns: Pi*sqrt(2)/2/GAMMA(3/4), which is about 1.8
pdf := z -> exp(-z^4)/c;
# now pdf integrates to 1
# we need the Statistics package for the next command,
# which was loaded above:
zDist := Distribution(PDF = pdf);

Now we can sample as follows:

# we will need the Statistics package, loaded above,
# and the distribution zDist:
xDist := Normal(0, sqrt(2)/2);
xValues := Sample(xDist, 100);
yValues := Sample(xDist, 100);
zValues := Sample(zDist, 100, method=[envelope,range=-2..2]);
XYZ := < xValues, yValues, zValues >^%T;

It's a little unfortunate that Maple can't decide by itself what the appropriate range for z is; the argument method=[envelope, range=-2..2] really shouldn't be necessary. I'll look into fixing that. I determined that -2 .. 2 would be a good enough range by evaluating:

CDF(zDist, -2.);
# returns about 1.9e-9.

Hope this helps!

Hi itsme,

The integration scheme used is Euler-Maruyama. This is implemented in QuantLib (we use version 0.3.11, by the way) in what their documentation calls the Monte Carlo framework -- in particular, in the PathGenerator<...> class. Their documentation page for this framework contains the following fragment:

The Black-Scholes equation

diff(f,t) + sigma^2/2*diff(f,x,x) + nu*diff(f,x) - r*f = 0

where r is the risk-free rate, sigma is the volatility of the underlying, and nu = r - sigma^2/2, has the form of a diffusion process. According to this heuristic interpretation, paths followed by the logarithm of the underlying would be Brownian random walks with a constant drift nu per unit time and a standard deviation sigma*sqrt(T) over a time T.

Therefore, the paths to be generated for a Monte Carlo model of the Black-Scholes equation will be vectors of successive variations of the logarithm of the underlying price over M consecutive time intervals Delta t_i, i = 0 ... M-1. Each such variation will be drawn from a Gaussian distribution with average nu Delta T_i and standard deviation sigma sqrt(Delta T_i) --- or possibly nu_i Delta T_i and sigma_i sqrt(Delta T_i) should nu and sigma vary in time.

That sounds like a pretty decent explanation of Euler-Maruyama to me.

Just to clarify: Wikipedia seems to suggest that Euler-Maruyama means that all time intervals are the same, but that's not the case for QuantLib: you can use any sequence of time points you like.

Hope this helps,

Erik Postma
Maplesoft.

Hi Bornin1992,

If you need to know the probability of more complicated events, you can use the ?Statistics/Probability command:

with(Statistics):
X := RandomVariable(Binomial(10, 1/2)):
Y := RandomVariable(Geometric(2/3)):
Probability(X + Y < 7);
# returns an exact fraction
Probability(X + Y < 7, 'numeric');
# returns 0.723...
# now to verify:
s := Sample(X + Y, 10^6):
SelectInRange(s, -infinity .. 6.999);
# should show a vector of roughly 723000 elements

Hope this helps,

Erik Postma
Maplesoft.

Hi GaoCG,

You'll need to provide initial values for your parameters so that Maple has a starting point for iterating towards a locally optimal solution. You can do this with the initialvalues option explained on the Statistics/Fit help page. Also, there seems to be an extra parenthesis in the call to Fit; I assume that's the one just before ', X, Y, sigma', right?

It's not so easy to find values of the parameters that lead to numbers that are anywhere near reasonable. Maybe you know the actual value of the parameters for some other instance of the model? Since there's a factor of sigma^2-(650*10^6)^2 in the denominator, where sigma^2 is on the order of a few hundred, this is a big negative number; since we need to find values that are positive, we need to compensate for that. One option would be to set A to a big negative number.

Another thing to remember is that there might be bad cancellation going on: because you're subtracting such a big number from sigma^2, if you want to see the effect of sigma on the result, you'll need sufficiently many digits to represent it accurately. It might be useful to increase Digits to a higher number, like 25.

Hope this helps,

Erik Postma
Maplesoft.

Hi hossayni,

What does currentdir(); print in these two worksheet? Maple will search for .m files in that directory, and if it doesn't match you either need to specify the path more fully, or change the current directory. See ?currentdir.

Hope this helps,

Erik Postma
Maplesoft.

Hi strawfields93,

The issue is in your use of 2n - if you use 2*n (in all three places where it occurs) it should work. (It did for me, at least.)

HTH,

Erik Postma
Maplesoft.

1 2 3 4 5 6 7 Last Page 1 of 10