Nash - Scientific Computing with PCs (523165), страница 46

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 46)

Several other packages were on ourshelves but not installed (they use up disk space!). We could also think of solving this problem in aspreadsheet package such as QUATTRO (our version of Lotus 123 does not offer multiple regression).Statistical packages offer extra output relating to the variability of the estimated parameters b as standarderrors and t-statistics, along with measures of fit and (optional) graphical checks on how well the modelworks. For purposes of illustration, we include a MINITAB output in Figure 15.5.5. Note that this wascreated by "executing" a file of commands prepared with a text editor.

Moreover, the output was savedby instructing MINITAB to save a console image (or log) file.The height difference problem is less straightforward. We can attempt to apply MINITAB, using aNOCONSTANT sub-command because X does not have a column of 1s. However, we do not know if thesoftware will cope with the indeterminacy in the heights (the "sea-level" issue). We can also use MATLAB’scapability to generate a generalized inverse of X. Data refinement problems like this one often arise inlarge scale land survey and geodesy calculations.

We know of a problem where X was 600,000 by 400,000.The structure of matrix X is very simple and the essential data for each height difference is just a rowindex, a column index, and an observed height difference. We do not need to store the matrix itself. If theproblem were large, a sparse matrix method is appropriate. A simple one is the conjugate gradient method(Nash J C, 1990d, Algorithm 24), though better methods could be found.15.5 Some SolutionsThis section presents some results mentioned above. For the weed growth problem, Figure 15.5.1 showsoutput from SYSTAT as well as the command file the user must supply to run this problem.

Figure 15.5.3gives the output of the Nash and Walker-Smith MRT code. Figure 15.5.2 is a graph of the data and thefitted model from the latter solution. We see from this, that the sigmoid shape of the logistic function isnot evident. This explains some of our difficulties. These are made clear in Figure 15.5.4, which gives theHessian and its eigenvalues at a point "near" the solution.A useful feature of the Marquardt (1963) method (Figure 15.5.3) is that it allows us to attach quite detailedpost-solution analysis to the program to compute approximations to standard errors and correlations ofthe estimated parameters. The user, in the problem specification, is also permitted in the Nash andWalker-Smith codes to include a results analysis that pertains to the application.

For example, we maywish to scale the parameters so the results are expressed in units familiar to other workers in the area.Figure 15.5.1Weed growth curve model using SYSTAT.a) Input command fileNONLINOUTPUT E:SSTTHOBBPRINT = LONGUSE HOBBSMODEL VAR = 100*U/(1+10*V*EXP(-0.1*W*T))ESTIMATE / START=2,5,3b) Screen output saved to a fileITERATION01. .

.8LOSSPARAMETER VALUES.1582D+03 .2000D+01 .5000D+01 .3000D+01.2889D+01 .2056D+01 .4985D+01 .3077D+01.2587D+01 .1962D+01 .4909D+01 .3136D+01DEPENDENT VARIABLE ISSOURCEREGRESSIONRESIDUALTOTALCORRECTEDVARSUM-OF-SQUARESDFMEAN-SQUARE24353.2132.58724355.7839205.4353912118117.7380.287RAW R-SQUARED (1-RESIDUAL/TOTAL) = 1.000CORRECTED R-SQUARED (1-RESIDUAL/CORRECTED) =PARAMETERUVWESTIMATE1.9624.9093.136A.S.E.0.1140.1700.069LOWER1.7054.5262.9801.000<95%> UPPER2.2195.2933.291ASYMPTOTIC CORRELATION MATRIX OF PARAMETERSUVWU1.0000.724-0.937V1.000-0.443W1.000134Copyright © 1984, 1994 J C & M M NashSCIENTIFIC COMPUTING WITH PCsNash Information Services Inc., 1975 Bel Air Drive, Ottawa, ON K2C 0X1 CanadaCopy for:Dr. Dobb’s JournalFigure 15.5.2Data and fitted line for the weed infestation problem.Figure 15.5.3Solution of the weed infestation problem by the Nash / Marquardt methodDRIVER -- GENERAL PARAMETER ESTIMATION DRIVER- 851017,85112821:52:4208-16-1992HOBBSJ.JSD 3-PARAMETER LOGISTIC FUNCTION SETUP- 851017FUNCTION FORM= 100*B(1)/(1+10*B(2)*EXP(-0.1*B(3)*I))HOBBSJ.JSD 3 parameter logistic fitbounds on parameters for problem0 <=b( 1 ) <=1000 <=b( 2 ) <=1000 <=b( 3 ) <=30ENTER INITIAL VALUES FOR PARAMETERS ([cr]=Y)B( 1 )= 1B( 2 )= 1B( 3 )= 1are masks or bounds to be set or altered?([cr] = no)MRT -- MARQUARDT/NASH -- 860412ITN 0EVAL’N 1 SS= 10685.29parameters111LAMDA = .00004LAMDA = .0004LAMDA = .004ITN 1EVAL’N 4 SS= 1959.699parameters.87460422.3246134.944148LAMDA = .0016LAMDA = .016ITN 2EVAL’N 6 SS= 1420.665parameters1.0035322.4409653.1544LAMDA = .0064...ITN 9EVAL’N 13 SS= 2.587254parameters1.9618614.9091623.135698LAMDA = 1.048576E-05LAMDA = 1.048576E-04LAMDA = 1.048576E-03LAMDA = 1.048576E-02ELAPSED SECS= 4 AFTER 10 GRAD & 16 FN EVALCALCULATED FUNCTION MINIMUM = 2.587254POST-SOLUTION ANALYSIS FOR MARQUARDT-NASH NLLSSIGMA^2 = LOSS FUNCTION PER DEGREE OF FREEDOM= .2874726PARAMETER CORRELATION ESTIMATESROW 1 :1.00000ROW 2 :0.72030 1.00000ROW 3 : -0.93651 -0.43731 1.00000SOLUTION WITH ERROR MEASURES AND GRADIENT OFSUM OF SQUARESB(1) = 1.961861 STD ERR = .1130655GRAD(1) = 6.923676E-04B(2) = 4.909162 STD ERR = .1688398GRAD(2) = -1.842976E-04B(3) = 3.135698 STD ERR = 6.863073E-02GRAD(3) = 8.602142E-04radii of curvature for surface along axialdirections& tilt angle in degreesfor B( 1 ) R.

OF CURV. = 1.580E-04tilt =-0.06543for B( 2 ) R. OF CURV. = 2.322E-03tilt =2.01501for B( 3 ) R. OF CURV. = 1.124E-04tilt = -17.24985RESIDUALS.011899 -3.275681E-02 9.203052E-02 .2087812.3926354 -5.759049E-02 -1.105722.7157898-.1076469 -.3483887 .6525879 -.287559515: DATA FITTING AND MODELLING135Figure 15.5.5 gives the MINITAB solution of the farm income problem.

The results appear satisfactory,though we should check them by plotting residuals against predicted values and against the explanatoryvariables. From a statistician’s perspective, we would want to plot the distribution of residuals andexamine some extra information in the output, for instance, the small magnitude of the t-ratio forpetroleum products, but we will not continue this discussion here.For the height difference problem, MINITAB (Figure 15.5.6) unfortunately does not apparently get refinedheight differences. It gives a warning message The difficulty is that the matrix X is now singular becausewe have no "absolute zero" of height defined, and MINITAB has acted to resolve the lack of uniquenessin the solution.

MATLAB is more revealing, as in Figure 15.5.7. Note that we have carried out thecalculations from a file SURVEY.M created with a text editor. MATLAB warns us that the X matrix (calledA here) is singular (that is, rank deficient). We regularize the MATLAB solution by subtracting b1 from allthe elements of b. This is equivalent to making benchmark 1 the "zero" of height.

Note that as part of our"solution" we compute residuals from the true heights of both our initial and final estimates, showing areduction in the sum of squares from 0.0016 to 0.0006. We also compute the errors from the true (butnormally unknown) heights used to create this problem. The error sum of squares is reduced from 0.0021to 0.00065. Generally we will not be able to know how well we have done.Figure 15.5.4THessian of the Hobbs weed growth nonlinear least squares problem for b = [200, 50, 0.3]a) Scaled parametersb) unscaled parametersInitial parameters?253Input scaling?100.010.00.1Hobbsf problem initiatedWARNING: sumsquares = 158.232at b253H indefinite :1.0e+004 *1.0314-0.25611.2851-0.25610.0645-0.34531.2851-0.34531.8578eigenvalues1.0e+004 *-0.00170.09492.8606End Warningeigenvalues of analytic Hessian1.0e+004 *-0.001701590222490.094854565720522.86059060485273Initial parameters?200.050.00.3Input scaling?111Hobbsf problem initiatedWARNING: sumsquares = 158.232at b 200.000050.0000 0.3000H indefinite :1.0e+006 *0.0000-0.00000.0013-0.00000.0000-0.00350.0013-0.00351.8578eigenvalues1.0e+006 *-0.00000.00001.8578End Warningeigenvalues of analytic Hessian1.0e+006 *-0.000000092166300.000000269649091.85779520216555Figure 15.5.5MINITAB output for the solution of the farm income problem (Table 15.1.1a)MTB > note REGDATA.XEC: regression of farm income index vsMTB > noteindex of phosphate and petroleum useMTB > read c1-c313 ROWS READMTB > end dataMTB > regress c3 on 2 c1 c2The regression equation isincome = 299 + 0.557 phosfate - 0.599 petrolPredictorConstantphosfatepetrolCoef299.00.55724-0.5989s = 10.37R-sq = 96.9%Analysis of VarianceSOURCEDFRegression2Error10Total12SOURCEphosfatepetrolStdev179.30.064220.8955DF11SS34135107535210SEQ SS3408748t-ratio1.678.68-0.67p0.1260.0000.519R-sq(adj) = 96.3%MS17067107F158.77p0.000136Copyright © 1984, 1994 J C & M M NashSCIENTIFIC COMPUTING WITH PCsNash Information Services Inc., 1975 Bel Air Drive, Ottawa, ON K2C 0X1 CanadaCopy for:Dr.

Характеристики

Список файлов книги