Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • My nl equation has an error, please help!

    Hello everyone,

    I am trying to estimate the parameters a and c of a nonlinear equation. The nonlinear equation is: y=a*(x/10)^(c-1). The code that I use is the following:


    Code:
    program nlces
    
        version 17
        
        syntax varlist(min=2 max=2) [if], at(name)
        local y: word 1 of `varlist'
        local x: word 2 of `varlist'
    
    
        // Retrieve parameters out of at
        tempname a c
        scalar `a' = `at'[1,1]
        scalar `c' = `at'[1,2]
    
    
        // Some temporary variables
        tempvar dterm
        generate double `dterm' = (`x'/10)^(`c'-1) `if'
        
    
        // Now fill in dependent variable
        replace `y' = `a'*`dterm' `if'
    
    end
    
     //Call nl like this, specifying initial values for some or all of the parameters:
    
    nl ces @ y x, parameters(a c) initial(a 0.0001 c 1)
    However, I get the error "nlces returned 9, verify that nlces is a function evaluator program". Could someone give me a hint on how to correct the error? In another question I read that the problem could be that I define the scalars as temporary, but that I refer to them in the -replace- function. Is this the problem? If yes, how can I solve it?

    Thank you very much for your help!

  • #2
    I would just go


    Code:
    nl (y = {a} * x^{c})
    and naturally you can apply extra arithmetic afterwards. Further, there is always a question of whether you would be better off taking logarithms to linearise that equation.

    Comment


    • #3
      Nick Cox Thank you very much! When I only apply your code, I get the error: "Iteration 0: cannot calculate derivatives equation/system not identified". Do you know why is this happening? Thank you in advance!

      Comment


      • #4
        I think in the equation that Nick has written the parameters a and c are not identified separately. For example if x=2, then both 2*2 = 4, that is a=2 and c=1; and 1*2^2=4, that is a=1 and c=2.

        Comment


        • #5
          And if I am wrong about the identification issue, you can give starting values like you did in your #1 as follows:

          Code:
           
           nl (y = {a=0.0001} * x^{c=1})

          Comment


          • #6
            I can't reproduce or even follow any of the supposed problems raised in #3 or #4 or #5, .


            Code:
            . clear
            
            . set obs 100
            Number of observations (_N) was 0, now 100.
            
            . set seed 2803
            
            . range x 1 100
            
            . gen y = 3 * x^0.7 + rnormal(0, 1)
            
            . nl (y = {a} * x^{c})
            
            Iteration 0:  residual SS =  39082.57
            Iteration 1:  residual SS =  28272.77
            Iteration 2:  residual SS =  23940.09
            Iteration 3:  residual SS =   20652.5
            Iteration 4:  residual SS =  15961.06
            Iteration 5:  residual SS =  13089.31
            Iteration 6:  residual SS =  128.2558
            Iteration 7:  residual SS =  102.3461
            Iteration 8:  residual SS =  94.46774
            Iteration 9:  residual SS =  94.46774
            
            
                  Source |      SS            df       MS
            -------------+----------------------------------    Number of obs =        100
                   Model |  237913.27          2  118956.634    R-squared     =     0.9996
                Residual |  94.467737         98  .963956498    Adj R-squared =     0.9996
            -------------+----------------------------------    Root MSE      =   .9818129
                   Total |  238007.74        100  2380.07735    Res. dev.     =   278.0965
            
            ------------------------------------------------------------------------------
                       y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
            -------------+----------------------------------------------------------------
                      /a |   3.016182   .0613011    49.20   0.000     2.894532    3.137831
                      /c |   .6979471   .0048235   144.70   0.000     .6883749    .7075192
            ------------------------------------------------------------------------------

            Comment


            • #7
              thank you for your comments, the last code works well! One last question, how can I condition the value of my parameter c, for example, to be bigger than 2? is this possible? Thank you again!

              Comment


              • #8
                If you have to force the fit, it's a very poor fit. Otherwise you need to reparameterise. some hints at https://www.stata.com/support/faqs/s...nts/index.html

                Comment


                • #9
                  Nick Cox got it! I will take a look. Thank you!

                  Comment


                  • #10
                    Here is the supposed problem raised in #4
                    Code:
                    . clear
                    
                    . set seed 666
                    
                    . set obs 10000
                    number of observations (_N) was 0, now 10,000
                    
                    . gen x = rnormal()
                    
                    . gen y = 1*x^2 + rnormal()
                    
                    . gen yy = 2*x^1 + rnormal()
                    
                    . nl (y = {a}*x^{c})
                    (obs = 10,000)
                    
                    Iteration 0:  cannot calculate derivatives
                    equation/system not identified
                    r(481);
                    
                    . nl (y = {a=2}*x^{c=1})
                    (obs = 10,000)
                    
                    Iteration 0:  cannot calculate derivatives
                    equation/system not identified
                    r(481);
                    
                    . nl (y = {a=1}*x^{c=2})
                    (obs = 10,000)
                    
                    Iteration 0:  cannot calculate derivatives
                    equation/system not identified
                    r(481);
                    
                    . nl (yy = {a}*x^{c})
                    (obs = 10,000)
                    
                    Iteration 0:  cannot calculate derivatives
                    equation/system not identified
                    r(481);
                    
                    . nl (yy = {a=2}*x^{c=1})
                    (obs = 10,000)
                    
                    Iteration 0:  cannot calculate derivatives
                    equation/system not identified
                    r(481);
                    
                    . nl (yy = {a=1}*x^{c=2})
                    (obs = 10,000)
                    
                    Iteration 0:  cannot calculate derivatives
                    equation/system not identified
                    r(481);
                    
                    .

                    Comment


                    • #11
                      I do not think you can enforce such a constraint (c>2) with the methods explained in the reference Nick provided https://www.stata.com/support/faqs/s...l-constraints/

                      I think the idea of the reference is to find a function which has a range coinciding with your constraint, e.g., if you want to enforce c>0 you use the exp(.) function. So you need a function which has a range c>2.

                      What you can do is to give initial values satisfying your constraint, and to hope that the solver will converge to a solution satisfying the constraint. If the constraint is not binding, the unconstrained and the constrained maxima will coincide.

                      Originally posted by Viridiana Corona View Post
                      thank you for your comments, the last code works well! One last question, how can I condition the value of my parameter c, for example, to be bigger than 2? is this possible? Thank you again!

                      Comment


                      • #12
                        As Joro Notes
                        it is possible that X has negative values. and if that is the case, the function itself cannot be numerically estimated.
                        THat is why it reports that odd problem.

                        Comment


                        • #13
                          Well, the report of #7 suggests that this does not bite with real data, although I cannot follow why #3 was issued following the same suggestion in #2.

                          Across several sciences fitting power functions (power laws, if you like) is standard and not especially problematic. in most cases.

                          Otherwise

                          Code:
                          nl (y = {a} * x^max(3, {c}))
                          is legal syntax and its results should be compared with those of an unconstrained fit.

                          Comment


                          • #14
                            It would be nice if what you are showing works, Nick. But I am not sure it does.

                            If we take your example where the model is apparently identified and the unconstraint solution is c=0.7, and I try to impose the constraint that c>1, look what happens:

                            Code:
                            . clear
                            
                            . set obs 100
                            number of observations (_N) was 0, now 100
                            
                            . . set seed 2803
                            
                            . . range x 1 100
                            
                            . . gen y = 3 * x^0.7 + rnormal(0, 1)
                            
                            . . nl (y = {a} * x^{c}), nolog
                            (obs = 100)
                            
                            
                                  Source |      SS            df       MS
                            -------------+----------------------------------    Number of obs =        100
                                   Model |  237913.27          2  118956.634    R-squared     =     0.9996
                                Residual |  94.467737         98  .963956498    Adj R-squared =     0.9996
                            -------------+----------------------------------    Root MSE      =   .9818129
                                   Total |  238007.74        100  2380.07735    Res. dev.     =   278.0965
                            
                            ------------------------------------------------------------------------------
                                       y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
                            -------------+----------------------------------------------------------------
                                      /a |   3.016182   .0613011    49.20   0.000     2.894532    3.137831
                                      /c |   .6979471   .0048235   144.70   0.000     .6883749    .7075192
                            ------------------------------------------------------------------------------
                            
                            . . nl (y = {a} * x^max(1,{c}))
                            (obs = 100)
                            
                            Iteration 0:  residual SS =  3052.151
                            Iteration 1:  residual SS =  3052.151
                            
                            
                                  Source |      SS            df       MS
                            -------------+----------------------------------    Number of obs =        100
                                   Model |  36030.414          0           .    R-squared     =     0.9219
                                Residual |  3052.1515         99  30.8298129    Adj R-squared =     0.9219
                            -------------+----------------------------------    Root MSE      =    5.55246
                                   Total |  39082.566         99  394.773392    Res. dev.     =   625.6309
                            
                            ------------------------------------------------------------------------------
                                       y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
                            -------------+----------------------------------------------------------------
                                      /a |   .8333161   .0095456    87.30   0.000     .8143756    .8522566
                                      /c |          0  (constrained)
                            ------------------------------------------------------------------------------
                              Parameter c taken as constant term in model & ANOVA table
                            So I do not know what Stata is doing here. First c should not be constrained, and second it should not be constrained at 0.

                            Originally posted by Nick Cox View Post
                            Well, the report of #7 suggests that this does not bite with real data, although I cannot follow why #3 was issued following the same suggestion in #2.

                            Across several sciences fitting power functions (power laws, if you like) is standard and not especially problematic. in most cases.

                            Otherwise

                            Code:
                            nl (y = {a} * x^max(3, {c}))
                            is legal syntax and its results should be compared with those of an unconstrained fit.

                            Comment


                            • #15
                              My own position as in #8 is that forcing a fit is not desirable. However, you can (try to) do it. My example is not to be taken literally If the real power is 0.7, then insisting it's at least 3 is a real strain on the fit.

                              Similarly,

                              Code:
                               nl (y = {a} * x^(3 + exp({c})))
                              is an example of how re-parameterisation ensures a power at least 3, contrary to the speculation in #11. Again, I don't recommend it.

                              Comment

                              Working...
                              X