/* DATE 4/24/91 */
/*     NONLINEAR REGRESSION
       Francesco Danuso, Istituto di Produzione Vegetale,
           Universita degli Studi di Udine, Udine, Italy.
       Joseph Hilbe, STB; English adaptation & revision
*/

capture program drop nonlin
program define nonlin
version 2.1
/* capture log close */
mac def scelta=1
while %scelta!=0  {
   !cls
   disp in wh _n(5) "                     NONLINEAR REGRESSION"
   disp _n(2) " Menu: " in wh _col(22) "0." in gr " Exit"
   disp in wh _col(22) "1."in gr" Model Fit"
   disp in wh _col(22) "2."in gr" Frequency distribution of residuals (Yo-Yp)"
   disp in wh _col(22) "3."in gr" Graph Yp vs. Yo"
   disp in wh _col(22) "4."in gr" Graph of residuals vs. Yp"
   disp in wh _col(22) "5."in gr" Graph of residuals vs. Xi"
   disp in wh _col(22) "6."in gr" Graph of Yp vs. Xi" _n(1)
   disp " Choice: " _req(scelta)
   !cls
   if %scelta!=0  {
      disp _n(1)
      disp " Cases (1=all): " _req(condi)
   }

   if %scelta==1  {
      disp _n(1) " Variable Y   : " _req(vary)
      disp " List of Xi   : " _req(varx)
      disp " Function     : " _req(funz)
      disp " # parameters : " _req(np)
      disp _n(1)
      /* azzera  */
      capture confirm var y0
         qui if _rc==0 {replace y0=.}
         qui if _rc!=0 {gen double y0=.}
      mac def j=1
      qui {
         while %j<=%np {
            capture confirm var a%j
            if _rc==0  {
               replace a%j=.
            }
            if _rc!=0  {
               gen double a%j=.
            }
            mac def j=%j+1
         }
         capture drop ycal zx
      }
      /* valin  */
      mac def j=1
      while %j<=%np  {
         disp "Initial value for b%j" _req(c%j)
         mac def j=%j+1
      }
     /* stipara - calculation  */
      mac def gira="c"
      mac def iter=0
      while "%gira"!="u" {
      if "%gira"=="p"  {
         /* valin */
         mac def j=1
         while %j<=%np  {
            disp "Initial value for b%j" _req(c%j)
           mac def j=%j+1
         }
      }
   mac def kk=%iter
   disp "How many iterations:" _req(nuit)
   disp _n(1)
   while %iter<%nuit+%kk  {
      /* trasb - transfer parameter values */
      mac def j=1
      while %j<=%np  {
         mac def h5="c%j"
         mac def b%j=%%h5
         mac def j=%j+1
      }
   /* calcin - parameter increasing for calc of partial deriv  */
      mac def j=1
      while %j<=%np  {
         mac def h6="c%j"
         mac def in%j=%%h6/10000
         mac def j=%j+1
      }
      quietly replace y0=%vary-(%funz) if %condi
      mac def i=1
      mac def varind=""
      while %i<=%np  {
         mac def h1="b%i"
         mac def h2="in%i"
         mac def b%i=%%h1+%%h2
         qui replace a%i=((%funz)-(%vary-y0))/%%h2 if %condi
         /* trasb  */
         mac def j=1
         while %j<=%np  {
            mac def h5="c%j"
            mac def b%j=%%h5
            mac def j=%j+1
         }
         mac def varind "%varind a%i"
         mac def i=%i+1
      }
/*   calculation of b  */
      quietly regress y0 %varind if %condi, noconst
      disp in gr " Iteration n.  " %iter+1
      mac def j=1
      while %j<=%np  {
         mac def h3="b%j"
         mac def h4="a%j"
         mac def b%j=%%h3+_coef[a%j]
         mac def c%j=%%h3
         mac def h10="c%j"
         disp " B%j= " %%h10
         mac def j=%j+1
      }
      mac def iter=%iter+1
  }
  disp _n(1)
  disp " Other initial values ->p; Out ->u; Continue ->c :" in wh _req(gira)
}
      /*  stampe  */
       disp _n(1)
  disp " Print of the results (y/n): "in wh _req(stampa)
  !cls
  if "%stampa"=="y" {
     log using prn:
  }
  quietly  {
    gen ycal=%funz if %condi
    disp in wh _n(1) "           =====NONLINEAR REGRESSION RESULTS====="
    disp _n(1)
    mac def iter=%iter
    disp " File: %S_FN" _col(45) "N. of iterations: %iter"
    disp " Variable Y  : %vary"
    disp " Variables Xi: %varx"
    disp " Model: %vary=%funz"
    disp " Data selection: if %condi" _n(1)
    summa y0, detail
    mac def resav=_result(3)
    mac def devst=sqrt(_result(4))
    mac def asim=_result(14)
    mac def curt=_result(15)
    disp "Residual Statistics:"
    disp "Residual Average = " %resav _col(35)  "Stand. Dev. = " %devst
    disp "Skewness         = " %asim  _col(35)  "Kurtosis    = " %curt
   regress %vary %varx if %condi
   mac def dtot=_result(2)+_result(4)
   mac def gltot=_result(3)+_result(5)
   mac def vtot=%dtot/%gltot
   disp _dup(60) "-"
   disp in gr "  Variation       d.f.         SS                 MS  "
   disp _dup(60) "-"
   regress ycal %varx if %condi, noconst
   mac def dmo=_result(2)+_result(4)
   regress y0 %varind if %condi, noconst
   mac def glmo=_result(3)
   mac def dre=_result(4)
   mac def glre=_result(5)
   disp in gr "  Model   " in ye _col(20) %glmo _col(30) %dmo _col(47)  /*
                            */  %dmo/%glmo
   disp in gr "  Residual" in ye _col(20) %glre _col(30) %dre _col(47)  /*
                            */  %dre/%glre
   disp in gr "  Total   " in ye _col(20) %glmo+%glre _col(30) %dmo+%dre /*
                            */  _col(47) (%dmo+%dre)/(%glmo+%glre)
   disp in gr "  Corr Total" in ye _col(20) %gltot _col(30) %dtot _col(47) /*
                            */  %vtot
   mac def rq=(%dtot-%dre)/%dtot
   mac def rq=int(%rq*10000+.05)/10000
   disp _dup(60) "-"
   disp in gr "  R^2 = " in ye %rq   _col(25)
   disp in gr _n(4)
   !cls
   mac def i=1
   /* standard errors  */
   while %i<=%np  {
      quietly test a%i
      mac def se%i=abs(_coef[a%i]/sqrt(_result(6)))
      mac def i=%i+1
   }
   disp _n(2)
   /* print coefficients */
   disp  _dup(70) "-"
   disp in gr "       Parameter        Standard Error           t           Prob. t"
   disp _dup(70) "-"
   mac def i=1
   while %i<=%np {
      test a%i
      mac def h1="c%i"
      mac def h2="se%i"
      mac def tst=%%h1/%%h2
      mac def h3=tprob(%glre,%tst)
      disp in gr " b%i  " in ye _col(13) %%h1 _col(30) %%h2 /*
                     */_col(47) %tst _col(62) %h3
      mac def b%i=int(%%h1*1000)/1000
      mac def i=%i+1
   }
   disp _dup(70) "-"
}
disp _n(2) " === CORRELATION COEFFICIENT AMONG PARAMETERS ==="
correlate, _coef
disp _n(6)
capture log close

      quietly gen zx=%vary if %condi
   }

   if %scelta==2  {
      !cls
      set textsize 130
      mac def salva ""
      mac def regi "n"
      disp _n(1)
      disp "Save graph (y/n)? " _req(regi)
      if "%regi"=="y" {
         disp "Graph name: " _req(salva)
         mac def salva "saving (%salva)"
      }
      mac def k1 "b1("        DEVIATIONS (Yo-Yp)") b2(" ")"
      mac def k2 "l1("Frequency")"
      gr y0 if %condi, hist bord xlab ylab bin(20) norm xline(0) /*
                   */ %k1 %k2 %salva
   }

   if %scelta==3  {
      !cls
      set textsize 130
      mac def salva ""
      mac def regi "n"
      disp _n(1)
      disp "Save graph (y/n)? " _req(regi)
      if "%regi"=="y" {
         disp "Graph name: " _req(salva)
         mac def salva "saving (%salva)"
      }
      mac def k1 "l1("PREDICTED VALUES")"
      mac def k2 "t1("%vary=%funz") t2("r^2=%rq")"
      mac def k3 "title("         MEASURED VALUES")"
    gr ycal zx %vary if %condi, bor s(Oi) c(.l) xla yla sort %k1 %k2 /*
                     */ %k3 %salva
   }

   if %scelta==4  {
      !cls
      set textsize 130
      mac def salva ""
      mac def regi "n"
      disp _n(1)
      disp "Save graph (y/n)? " _req(regi)
      if "%regi"=="y" {
         disp "Graph name: " _req(salva)
         mac def salva "saving (%salva)"
      }
      mac def k1 "b1("    PREDICTED VALUES")"
      mac def k2 "t1("%vary*%funz") t2("r^2=%rq")"
      mac def k3 "l1("    DEVIATIONS (Yo-Yp)")"
      gr y0 ycal if %condi, bord s([_n]) xla yla yline(0) %k1 %k2 %k3 %salva
   }

   if %scelta==5  {
      !cls
      set textsize 130
      mac def salva ""
      mac def regi "n"
      disp _n(1)
      disp "Save graph (y/n)? " _req(regi)
      if "%regi"=="y" {
         disp "Graph name: " _req(salva)
         mac def salva "saving (%salva)"
      }
      describe %varx
      disp "Variable Xi: " _req(vx)
      mac def k1 "l1("    DEVIATIONS (Yo-Yp)")"
      gr y0 %vx, bord s([_n]) yline(0) xlab ylab psize(70) %k1 %salva
   }

   if %scelta==6  {
      !cls
      set textsize 130
      mac def salva ""
      mac def regi "n"
      disp _n(1)
      disp "Save graph (y/n)? " _req(regi)
      if "%regi"=="y" {
         disp "Graph name: " _req(salva)
         mac def salva "saving (%salva)"
      }
      mac def k5 "t1("Symbol diameters proportional to (hidden) %vna")"
      describe %varx
      disp "Variable Xi plotted: " _req(vx)
      disp "Variable Xi hidden : " _req(vna)
      quietly summ %vna if %condi
      quietly gen peso=99*(%vna-_result(5))/(_result(6)-_result(5))+1
      gr %vary %vx =peso, bord s(O) xlab ylab psize(60) %k5 %salva
      drop peso
   }
}
set textsize 120
end