In the case of multiple solutions to solve(diff(func, x)=slp, x) you can modify the "g" function in my previous code by this one

g := proc(func, slp)
  local a, b;
  a := select(type, [solve(diff(func, x)=slp, x)], 'realcons');
  b := [seq(y[k] = Tangent(func ,x=a[k]), k=1..numelems(a))]:
  print(func, slp, 1);
  Maplets:-Tools:-Set('MMLV1'=b)
end proc:

Tip: The select(...) command is necessary because Tangent requires a[k] to be real ; look to what happens if func=x^4 if the line a:=... is replaced by a := solve(diff(func, x)=slp, x);

The output in the MMLV1 field is of the form [tangent1, tangent2, ...] and not very convenient (my opinion): writing each solution on different lines would probably looks better (?).
Unfortunately MathMLViewer doesn't enable this (other solutions do exist)

@epostma 

I use Maple 2015

But, damn it! I realized my test with Student:-Statistics by writing  
R := NegativeBinomialRandomVariable(2, 0.003);
instead of 
R := NegativeBinomialRandomVariable(1.965, 0.003).

Here is, realized with Maple 2015, a quick comparison between the two packages.

NegBin_2015.mw

They both behave de same way (giving identical results on the same test cases or failing simultanously).
Then  my sentence "I discovered inadvertently that the package Student[Statistics] worked perfectly on the same problem, even for a success probability as small as 3e-10" is completely wrong.
Apologies.

My writing  R := NegativeBinomialRandomVariable(2, 0.003);  comes from the fact that a number of failures of 1.965 has no meanings: this number must me a positive integer.
The same occurs with the Binomial(n, p) distribution which accepts on integer values for n.
To be honest Maple is not the only language to admit this largesse which, from a probabilistic point of view is a nonsense.
This probably comes from the inner representation of Binomial/NegativeBinomial in terms of the Gamma function (defined over R) instead of true factorials (even if the Gamma function can be thought of an extension of factorials, they are not the same thing).
 

for info, 

p := 0.05;
alpha:= 0.98;
Digits := 20:
_EnvStatisticsIterations := 10000:
Q := Quantile(X, alpha, numeric);
Probability(X > Q);

                              0.05
                              0.98
Warning, numeric quantile computation failed to converge; consider increasing _EnvStatisticsIterations or Digits.
                              FAIL

The default value of _EnvStatisticsIteration is 100: even increasing by a factor 100 anfd pushing Digits uo to 20 doesn't solve the problem.


Q := Quantile(X, alpha, inert);
value(%);

fails too.

Finally: use package Student[Statistics] instead of package Statistics
restart:
with(Student[Statistics]):
X := NegativeBinomialRandomVariable(2, 0.003):
Quantile(X,0.98, numeric);

                             1941.
X := NegativeBinomialRandomVariable(2, 3e-10):
Quantile(X,0.98, numeric);
                  1.944640566866634e10

@tomleslie 

Very rude reply, a simple copy paste can do the job with minor other actions.
I don't think it was necessary to answer this way.

I you have constructive ideas to propose, here is the code you requested


 

Download NegBin.mwNegBin.mw

Tha package Statistics has very interesting features and at the same time suffers a large number of weaknesses that make it unusefull for even very simple things
In particular Maple doesn't manipulate correctly discrete random variables in many situations.
You just faced one.

To show you what a software dedicated to statistical analysis does, here is a screen capture obtained with the language R

R.pdf




The main difference between Maple and a software dedicated to statistical analysis, is that 

@acer 

Thanks acer, 

It's just a part of a more complicated problem it would take too much time to detail.

In a few words:
Let F an application from the unit cube K into R²,  F(K) the image of K by F and  B a subset of F(K).
My goal is to find the subset A in K such that F(A)=B

F is extremely complicated (not even known in anlytical form).
To handle the problem I developped a procedure wich returns a set of points X = {x[n], n=1..N} which verify F(x[n]) in B for each n=1..N.
This procedure starts from a set points randomly choosen in K and evolves to the set X. I assume that A can be reasonably estimated by the convex hull of X
Now I need to "visualize" A to understand the structure of F.
My first attempt was to use simplex:-convexhull to visualize the (approximated) projections of A onto the three canonical plans of the input space.

Next I discovered the PolyhedralSets package that could enable to represent the 3D convex hull of X.
But, as I said, PolyhedralSet (first step before using ConvexHull) builds the polyhedral which vertices {x[n], n=1..N} ONLY if the coordinates ox eaxh x[n] are rational.

So my question.

Your suggestion about ComputationalGeometry:-ConvexHull is very interesting (I wasn't aware of this package).
I will try it in a few days as soon as I will return to my office.

Than you acer

@Rohith 

Use the astute vv's idea in a slightly different form


 

Download Inspired_by_vv.mw

@tomleslie 

TimeSeries corrected

TimeSeries.mw

I was right, could have been done by redefining tickmarks... but it's better to do the job correctly.
Then, contrary to what you wrote "... your reference to your earlier  TimeSeries.mw is invalidated by the simple fact that the OP's data is in 2-year steps ..." is too crude an affirmation


More of this, if I wanted to be fussy I would point ou to you that if x[1] is 1978, then x[30] is not 2038 but 2036.

The devil is in the detail and nobody' perfect 

@tomleslie 

"The simple reason I don't propose to add much further to this thread is that I am not a great believer in lengthy extrapolation from actual data. Such extrapolation is always mathematically possible - but tiny discrepancies in fitted functions can lead to huge discrepancies in extrapolated values, so is the process meaningful/worth it"

I agree one hundred 100%!

I haven't payed attention to the "count by two's".
I have no time right now to correct my analyses, but unless I'm mistaken, it will just correspond to a change in the horizontal abscissa. Just if I wrote tickmarks=[1978=1978, 1979=1980, 1980=1982, ...]
I will redefine the time serie by writing explicitely 'dates' = ["January:1978:01", "January:1979:01", ...]  to be sure.

Hope to see you again on a future thread

It's not possible to compute the coefficient of sinh(n*Pi/lambda)  in the expression g
Simply because g is not a polynomial of the indeterminate  sinh(n*Pi/lambda)  (first line of the coeff help page)


 

Download coeff.mw

@tomleslie 

Probably my last contribution (all good things must come to an end)
You will find in this file an estimation of the  prediction error of x[30] for the original x[n+1]=a*x[n]+b

Using the least squares method to estimate (a, b) through Statistics:-LinearFit returns the covariance matrix of the estimators of (a, b).
This covariance matrix is used to assessed the distribution of x[30].

NOTA: The classical formula ued to assess the prediction error of y[n] = a*x[n]+b  (note y, not x[n+1])
var(a*X+b)=MSE*(1/N+(X-(Mean(x))^2 / Variance(x) ) is of no use here.

Variability_of_the_prediction.mw

@tomleslie 

The underlying model constructed in my file (with ExponentialSmoothingModel) is not the simple model the OP suggested.
I know from my student years what exponential smoothing means but its implementation in Maple is far more complex than what I was taught years ago.

ExponentialSmoothingModel has a lot of options, that's why I used before the command Specialize plus AIC, BIC, Optimize, in order that Maple returns the "best" specialized exponential smoothing model.
Each of these three criteria say that  < an ETS(M,M,N) model > is the best model, that is:

• the evidence stands for multiplicative errors
• and multiplicative trend too
• there is no evidence for a seasonal effect
   

Note that we could force the model to account for additive errors.

One the "best" (better to say "the more likely") model identified, I used OneStepForecasts to forecast x[n+1] from x[n] (first plot).
This plot is completely different from the plots you provided (and I did too in my earlier replies) for the "model" curve is not a straight line but a part of an exponential curve (or kind of).
This is the major reason why the forecast of x[30] is higher (mean value 3.94) than the one obtained by the x[n+1]=a*x[n]+b original model (second figure).

To obtain x[30] (more precisely x[k], k=14..30) I used Forecast. This function enables computing confidence intervals too.
The 95% bilateral confidence interval for x[3] is about 3.94 +/- 1, and thus "rather" large.
We could expect that becasu we "extrapolate" very far from the last point x[13], because the first plot shows rather large residuals errors, and because these later are amplified by the exponential model.

In a nutshell, even if this exponential smoothing model is not the one proposed by the OP, I would prefer the former from a prediction perspective (more solid basis) ... until I had strong prior reasons to consider that x[n+1]=a*x[n]+b is "the good model".
If I was faced to this time serie, without prior knowledge about a reliable underlying model, I would use this TimeSeries package, with the minimum bunch of precautions (use of Specialize with AIC or BIC).
Of course, if someone has very small knowledge about time series ans the way to handle them, TimeSeries will probably look like an unacceptable black box.


That is the meaning of my " ...a package such as TimeSeries is interesting: you can do (and I can too) do more reliable things by using it"

@tomleslie 

Hi,

I think that phe package TimeSeries would be more appropriate, at least to answer the Zeineb's s second question about the forceasting of x[30].
I sent a reply to Zeineb some hours ago. I'm not very familiar with the TimeSeries package, so I think improvements can be done.
You wrote: "the isue with auto regressive ...". Surely, that's why a package such as TimeSeries is interesting: you can do (and I can too) do more reliable things by using it


----------------------------------------------------------------------------------
For info: how can we choose among a linear or a quadratic model?
A very simple solution is based upon the so-called "generalization error", by opposition to the classical "representation error" (plus or minus the residual sum of squares in classical linear regression)
Obviously the higher the degree of the fitting polynomial, he lower the representation error (null when the degree equals the number of observations minus 1).
But as you know high degree polynomials use to oscillate strongly in some situation.
A polynomial that strongly oscillates between observed points is said to have a large generalization error (that is a potentially large error far from the observed points).
The idea is that, between the 0 order degree polynomial (maximal representation error, no oscillations)) and the 12 degree one (minimal representation error, and potentially large generalization error), there exist a compromise between a not to large representation error and a rather small generalization error.
To asses the latter one uses to use "resampling methods".
All of them is based on the following ides: one uses only a part of the original data to fit a model of a given order, and we asses its generalization error by running it on the remaining part of the data.
The simplest method is Leave-One-Out, look at  https://en.wikipedia.org/wiki/Leave-one-out_error  for instance.
A more sophisticated approach is the Bootstrap method that the Statistics package already proposes.
The "better compromise" is generally found by computing the AIC or the BIC information criterion.

Not your original request but this could nevertheless interest you.
Maple posseses a "TimeSeries" package to treat the kind of forecast pronlem you have posted.

TimeSeries.mw

# prediction of x[30] from differents starting points from x[1] to x[12]

for k from 1 to N do
printf("%2d : %1.5f\n", k, (x@@(30-k))(data[k]));
end do;
 1 : 1.56716
 2 : 1.55764
 3 : 1.41635
 4 : 1.40748
 5 : 1.39675
 6 : 1.55937
 7 : 1.37024
 8 : 1.30034
 9 : 1.44284
10 : 1.47176
11 : 1.49616
12 : 1.61075
13 : 1.53250