mmcdara

6149 Reputation

17 Badges

9 years, 72 days

MaplePrimes Activity


These are replies submitted by mmcdara

@ogunmiloro 

Here are a few linear fits of these new data.

The problem is interesting for it enables introducing a basic trick in statistical regression: the scaling of the data.
The worksheet must be read step by step:

  1. In a first attempt I try to fit (Statistics:-LinearFit) a high degree (14) polynomial without any precaution.
    The result is obviously non acceptable (you can reduce the value of P [degree of the polynomial] to see when the "problem" begins to occur).
     
  2. The reason comes from the huge values the leading monomial being xP takes when X is close to 300. 
    This is a classical issue which can be (generally) fixed by centering and scaling the data.
    Once done Statistics:-LinearFit returns an acceptable model.
     
  3. The exact reason of the problem is detailed in step 3: the condition numbers of the "Fisher matrix" are computed for the RAW data (about 1070) and SCALED data (about 1011).
    The algebra used in Statistics:-LinearFit is very simple and I wrote the relations that give the values of the regression parameters (vector beta).
    Note I still operate with SCALED data.
     
  4. I show now how to transform the model obtained on SCALED data into a model on RAW data.
    From the plots we can make a few observations:
    1. the model presents some oscillations near X=0 and is not as flat as Y is,
    2. the model decreases near X=300 (this comes from the last value of Y which "attracts" the model (this is something thar cn be investigated with Statistics:-LinearFit [on SCALED data] by studying the leverages abd the Cook's distance [ output=[leverages, CookDstatistic] ]
       
  5. In this final step I show how to constraint the model in order it be flatter near X=0 and increasing for the higher values of X.
    The solution relies upon  Optimization:-Minimize where an objective function J (residual sum of squares) is minimized under constraints
    1. the first 5 derivatives of the model are null at XC[1] (which corresponds to X=0)
    2. the value of the model at XC[1] = YC[1]
    3. the first and second derivatives of the model are non negative on the las 10 points


Fit_of_your_NEW_data.mw



 

 

@ogunmiloro 

On point (2) the answer is "technically yes" but, as Carl noted it seems less "natural"... unless you have evidence to say that the expectation of the model is indeed a "double sigmoid"  and the data equal this trend plus an additive "noise" (which is the framework of statistical regression).
 

@Carl Love @ogunmiloro

Thanks Carl  for your comments
 

"Show the extrapolation of the sum-of-sigmoids fit out to its upper horizontal asymptote. I think x=500 is a good stopping point."
That is a good idea, so I have completed my previous reply accordingly.

One of my prefered type of model is named  Kriging.
It is implemented in Maple since, I believe, Maple 2019 to answer problems related to spatial statistics (aka geostatistics).
Thet are very versatile models, even extremely parcimonous (fhey may give astonishingly good models [if you use them carefully] even if they depend only on a small handull of parameters).

There are different forms of Kriging (ordinary, simple, universal, ...).
All of them are of the form a Trend plus a Gaussian Process
In the Universal Kriging (UK) the trend is generally a function that already captures long sacle behaviours, leaving it to the gaussian process component to capture small scales.
(In simple Kriging the trend is a constant and, in some sense, the gaussian process does almost all the job).

The idea is to use the 6 parameters (positive) model and to complete, or correct) it by a gaussian process
The final model is

  • positive for all x
  • monotonically increasing for almost all x
    (oscillations at earlier values of x comes from the fact that Y has  identival values repeated several times)
  • equal to 1 for X=1
  • bounded for extrapolation for x >> max(X)=298

The computation of the error of prediction is not given (the procedure KG is not a full implementation of a Kriging emulator).

On the left most figure zoom in the region defined by the rectangle to see what happens when extrapolation begins;
The Kriging model abruptly connects to the trend over a x-range which is of the order of the correlation length (chosen to 2.5).
 





Kriging_of_your_data.mw

@ogunmiloro 

I think that you would have said clearly what you wanted in your initial question.

What do you want:

  • statistical estimation of parameters of an a priori model
    • your a priori model A*exp(r*x) 
    • a polynomial model (tomleslie) which won't respect positivity constraints (that you seem to say now they are important!!!)
  • smoothing
    • spline interpolation model (Kitonum) (whose residual sum of squares is obviously 0)

Why do you want to obtain such a "fit"

  • to predict the future ("what Y can I expect for X > 298?"
    in this case spline interpolation or polynomial fitting is useless
  • to interpolate in the range X=0..298
    then use spline interpolatoin


Please look to my last reply.


 

@tomleslie 

This is not an answer to the question!

@Carl Love 

Thanks, excellent!

@Carl Love 

Nice way.
I vote up even if this doesn't work in Maple 2015 :-)
(Fr is correct)



It is a shame that DrawGraph transforms rationals into strings whereas TEXT can display rationals  (I followed the track with the library assistant without finding exactly where it's done).

More of this I discovered that only rational or float weights are enabled by Graph (weights like Pi, sqrt(2), ... generate an error).

@Carl Love 

There's a silver lining to everything: writing snippets of code (I do not remember there did exist a few years ago) probably leads to smaller contributions (some take almost a whole mapleprimes' page), attaching an image can be more explicit and appealing than putting it in the middle of a long worksheet ... 
Personally I think it's not that stupid to favour code snippets instead of  "file contents", but if it is a policy I feel it must be stated clearly.

Wait and see

@dnaviaux 

No problem

Hi, 
Could you detail what you mean by "table":

  • is it a "table" in the Maple's sense
  • or the type of "table" one get using DocumentTools:-Layout for instance

Is this question related to this thread What might be the best way to convey calculated re... ?

 

@dnaviaux 

By the way, I didn(t pay a careful attention to all the contributions about the "presentation" of the results.
Maybe this alternative to your procedure rprintf might interest you

SEA_SI_3.mw

@dnaviaux 

Thnak you,
I will look at it with interest, ... after Christmas Eve.

Happy Christmas

 

@dnaviaux 

To be a little bit more exhaustive, here is a new worksheet which

  1. contains the previous one 
  2. gives an answer to calculate a kind of worst case (extremal values)
    1. first by using the Tolerances package... but as you already observed yourself it fails when the expressions are too complicated
    2. secondly by Monte-Carlo simulation

The coding I use for Monte-Carlo simulation could be enhanced when the quantity of interest is not a function of an independant variable (t in my toy example)  but, as it seems to be your case, the solution of an algebraic system.

SEA_SI_2.mw

@Carl Love 

You were right, the good term to use was "parameter", "factor" or even "contributing factor" to be more precise, but not "contributor" whose, as you said, "usually refers to a person who contributes".

Here again there is a subtle difference between correcting and editing.
As you can see here "to correct" has few translations in French.
https://www.linguee.fr/francais-anglais/search?source=auto&query=to+correct
but one of them ("corriger") has many more in English (including "to correct").
https://www.linguee.fr/francais-anglais/search?source=auto&query=corriger
In French, "correct" (= "to correct") doesn't make as many different meanings as in English, and "to correct" can be either nice or rough depending on the words you use with it.

 

 

 

@Carl Love @dnaviaux

 

Thanks Carl for having corrected me. 
As you've already noticed, I sometimes use the wrong term, and in this cas I'm all the more shameful that I lnew it perfectly well  "importance factor" or, as you wrote "factor that contributes to the error".
That's what happens if you think in your vernacular when writing in an other language.

In any case, here is @dnaviaux   a toy example which whows of to assess the effect of each factor.
Remember I use in a prcious reply the notional expression y=f(a, b, ...).
When all the (random variables which model) a, b, .. are mutually independent, the propagated variance writes
Variance(y) = diff(y, a)^2 * Variance(a) + diff(y, b)^2 * Variance(b) + ...

A classical definition of the "effect of a factor" (let's take "a" for instance) is 
Effect(a) = diff(y, a)^2 * Variance(a) / Variance(y) .

As SEA returns Variance(y), the method to obtain each "Effect" is simply to run SEA as many times as the number of factors while keeping fixed all the factors but the one of interest.
For instance:

  • ​​​​I write for short Variance(y) = SEA(a, b, c, ...)
  • Each a, b, c is defined this way : a := Quantity(ma, sa)), b := Quantity(mb, sb)), ...
  • Let a0 := Quantity(ma, 0)), b0 := Quantity(mb, 0)), ...
  • To assess the effect of factor a you just have to compute 
    SEA(a, b0, c0, ...) / SEA(a, b, c, ...)
     

SEA_SI.mw

 

PS 1 : I have said many times that the factors a, b, ... are modeled by random variables (when you use SEA for instance). This raise many questions:

  • If your Uncertainty is a constant value (lets say a = ma +/- Ua) then how are you going to define the value of sa when you define a as a:=Quantity(ma , sa)?
    • if you think the uncertainties are gaussian you will probably write something like 
      sa = Ua/2 or sa = Ua/3
    • if you believe they are uniform, then you will use sa = Ua/sqrt(3)
       
  • If your Uncertainty is a percentage (lets say a = ma *(1 +/- Ua)) this means the parameter  sin a:=Quantity(ma , sa) is not constant, which can be very complicated to handle if its mean (nominal) value ma varies in large proportions

PS 2 : If not all the funcertainties are mutually independent the formula for computing the effect of a factor is more complicated, but the more important point is that the definition of this effect is less simple.
The formula I gave above only accounts of what is called the "elementary effect" of (there) factor a.
If the uncertainties of a and b are correlated, then Variance(y) contains a term
2*diff(y, a)*diff(y, b)*Covariance(a, b).
This means that a contributes to Variance(y) through its interaction with b "interaction effect".
For a given factor there exist an "elementary effect" and a "total effect" wich includes the "elementary effect" plus all the effects through interactions this factor has with the other ones.

First 63 64 65 66 67 68 69 Last Page 65 of 125