2
$\begingroup$

I think that this is an interesting issue for all those who when to not only adjust one prediction function, but also to obtain the prediction bands to asses the graphical or tabular error I'm trying to apply step by step the Delta Method to my data, however, as I would like to share with you my advances to get to the final. Summarizing and reviewing the standard method using in commercial software above all, we find: Delta Method in confidence and prediction intervals

Where the right square root involves the inverse of the Hessian matrix I think, although they call it design matrix. As well as the gradient matrix function. So I share here my code to try to fit a custom non-linear function and assessing its prediction intervals (confidence and prediction, or int Mathematica language, mean and simple prediction intervals).

CODE WITH MY EXAMPLE DATA:

(*Example data*)
data = {31, 46, 70, 87, 87, 93, 114, 128, 133, 134, 143, 155, 161, 
161, 163, 177, 181, 207, 207, 226, 302, 315, 319, 347, 347, 362, 
375, 377, 413, 440, 447, 461, 464, 511, 524, 556, 800, 860, 880, 
954, 5200, 12000};

dist = ProbabilityDistribution[{"CDF", 
Exp[-a Exp[-x b] - c Exp[-x d]]}, {x, 0, Infinity}, 
Assumptions -> {{0 < a < 10}, {0 < b < 1}, {0 < c < 10}, {0 < d < 
1}}];

res = FindDistributionParameters[data, 
     dist, {{a, 3.8}, {b, 0.006}, {c, 0.08}, {d, 0.0002}}];

            (* ACOV is:*)
acov[data_, dist_, paramlist_, mleRule_] := 
     Block[{len, infmat, cov}, len = Length[data];
     infmat = -D[LogLikelihood[dist, data], {paramlist, 2}]/len /. 
     mleRule;
cov = Inverse[infmat]];

acov[data, dist, {a, b, c, d}, res];

           (* Gradient matrix G', in this case with only one variable x, we obtain a vector with only one element *)
In[1]:= D[Exp[-a Exp[-x b] - c Exp[-x d]], {x}]/. res;


            (* To multiply Transpose-Gradient by acov by Gradient, is not necesary to transpose the first matrix, Mathematica detect the same matrix and does it for us *)

%2.%1.%2

            (*T-student critical value to alpha 0.05 and 42 data*)
Quantile[StudentTDistribution[41], 0.95]

The question/s here are:

  1. What snips or lines of code I must to add to complete and apply my code to a numerical example as we can find in here, getting the predicted y and the their confidence and prediction intervals? As if you try my code it does't give the table of confidence intervals as does exp["MeanPredictionConfidenceIntervalTable"] in the case of Non-Linear fiiting (that is different of fitting a distribution with the MLE method). That is because I have some doubts and I am not an old user of Mathematica:

    1.1.: A little doubt is if finally my ACOV matrix is well assessed with my code. I don't understand why in the post I got the first information about the Hessian and Information matrices they divide the matrix of partial derivates by N (Length[data]). Please, check it.

    1.2.: Following the Delta method I have to execute the dot product of the Transpose-Gradient by ACOV by Gradient, that is, %2.%1.%2, so I obtain a simple number. Is this correct?. But I found that if we haven't the same number of columns in the vector that the number of rows in the matrix we obtain other matrix, not a number. If we fill those columns with 0, then we obtain a number, but maybe this isn't the correct way.

Please, comment this post if you have any is not clear or edit my code if you know the appropiate code lines.

Reference posts:

I look forward to your help, thank you.

$\endgroup$
  • 2
    $\begingroup$ Sadly, as it is currently written this is not a Q :) If this is intended to be a Q add any relevant details/code and try to keep it simple. $\endgroup$ – Sektor Jul 28 '14 at 10:18
  • 1
    $\begingroup$ I too find this post confusing. On the one hand it seems like you are setting up for a self-answered question to help others ("I would like to share with you my advances") yet no answer is given. On the other you are seeking help ("I look forward to your help") yet no clear question is asked. What is your intention? $\endgroup$ – Mr.Wizard Jul 28 '14 at 10:40
  • $\begingroup$ My appologies, I believed I could upload a file, and then I forgot to add my code. Now I added my snip of code. $\endgroup$ – JosGranada Jul 28 '14 at 11:22
  • $\begingroup$ Even with the code, it's not clear what the question is. But please do edit the question, because there might be an interesting question underneath it. $\endgroup$ – Verbeia Jul 28 '14 at 11:31
  • 1
    $\begingroup$ @bobthechemist If so it's from years of practice repeatedly correcting myself. :o) $\endgroup$ – Mr.Wizard Jul 28 '14 at 11:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.