(newcommand{xb}{{bf x}}
newcommand{betab}{boldsymbol{beta}})Earlier than you utilize or distribute your estimation command, you must confirm that it produces right outcomes and write a do-file that certifies that it does so. I focus on the processes of verifying and certifying an estimation command, and I current some methods for writing a do-file that certifies mypoisson5, which I mentioned in earlier posts.
That is the twenty-fifth publish within the collection Programming an estimation command in Stata. I like to recommend that you simply begin initially. See Programming an estimation command in Stata: A map to posted entries for a map to all of the posts on this collection.
Verification versus certification
Verification is the method of building {that a} command produces the right outcomes. Verfication produces true values that may be in contrast with the values produced by a command. Certification is the method of checking that the variations between the verified true outcomes and the outcomes produced by a command are small enough.
Verification might be simple or tough. If there may be one other command or program that you simply belief, you should use it to create verified values. For instance, I belief the poisson command, so I can use it to create true test-case values for mypoisson5. When one other command or program is just not out there, I exploit simulation methods to acquire licensed values. See Monte Carlo simulations utilizing Stata and Effectivity comparisons by Monte Carlo simulation for discussions of Monte Carlo simulations.
I certify a command by writing do-files that test that the outcomes of a command are near the verified true values in lots of particular circumstances. These do-files are referred to as certification scripts, and I run them each time I make any change to my command. The method that I current is a tremendously simplified model of that used to certify Stata; see Gould (2001) for one more introduction to certification and for extra about Stata certification.
Evaluating numbers
The assert command checks {that a} logical expression is true. Right here I exploit it to test that two integer values are equal or that two noninteger values are sufficiently shut.
I test for equality between integer values and closeness between noninteger values due to how computer systems do math. You can’t match your entire real-number line on a pc; there are too many actual numbers. Computer systems use finite-precision base-two approximations to the true numbers. Integers have a precise illustration on this approximation, and integer calculations might be carried out with out approximation error. Most noninteger values should not have a precise illustration within the base-two approximation used on computer systems, and noninteger calculations are carried out with approximation error. See The Penultimate Information to Precision for particulars.
Instance 1 illustrates that assert produces no output or error when requested to say {that a} true logical expression is true.
Instance 1: assert a real expression
. assert 3==3
In distinction, instance 2 illustrates that assert produces an error when requested to say {that a} false logical expression is true.
Instance 2: assert a false expression
. assert 3==2
assertion is fake
r(9);
In instance 3, I exploit mreldif() to compute the utmost of the element-wise relative variations between two integer-valued vectors, and I then use assert to test for equality.
Instance 3: assert and mreldif()
. matrix a = (1, 2, 3)
. matrix b = a
. show mreldif(a, b)
0
. assert mreldif(a,b)==0
Examples 1–3 illustrated assert by evaluating integers. In certification, we often evaluate noninteger values. Due to finite-precision approximation errors, a small change in how the outcomes are computed—akin to utilizing 1 processor as an alternative of 8 processors in Stata/MP, altering the kind order of the info, or utilizing a Mac as an alternative of a PC—might trigger the outcomes to vary barely. These adjustments might be shocking in case your instinct is guided by infinite-precision math, however they need to be sufficiently small to be ignorable. Instance 4 illustrates this level by evaluating the purpose estimates obtained utilizing 8 processors and 1 processor.
Instance 4: The impact of utilizing 1 as an alternative of 8 processors
. clear all
. use accident3
. gsort - cvalue
. set processors 8
The utmost variety of processors or cores getting used is modified from 1 to
8. It may be set to any quantity between 1 and eight
. mypoisson5 accidents cvalue children visitors
Iteration 0: f(p) = -851.18669
Iteration 1: f(p) = -556.66855
Iteration 2: f(p) = -555.81731
Iteration 3: f(p) = -555.81538
Iteration 4: f(p) = -555.81538
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .0706484 -9.28 0.000 -.7943553 -.5174188
children | -1.009017 .0807961 -12.49 0.000 -1.167374 -.8506596
visitors | .1467115 .0313762 4.68 0.000 .0852153 .2082076
_cons | .5743541 .2839515 2.02 0.043 .0178194 1.130889
------------------------------------------------------------------------------
. matrix b1 = e(b)
. set processors 1
The utmost variety of processors or cores getting used is modified from 8 to
1. It may be set to any quantity between 1 and eight
. mypoisson5 accidents cvalue children visitors
Iteration 0: f(p) = -851.18669
Iteration 1: f(p) = -556.66855
Iteration 2: f(p) = -555.81731
Iteration 3: f(p) = -555.81538
Iteration 4: f(p) = -555.81538
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .0706484 -9.28 0.000 -.7943553 -.5174188
children | -1.009017 .0807961 -12.49 0.000 -1.167374 -.8506596
visitors | .1467115 .0313762 4.68 0.000 .0852153 .2082076
_cons | .5743541 .2839515 2.02 0.043 .0178194 1.130889
------------------------------------------------------------------------------
. matrix b2 = e(b)
. show mreldif(b1, b2)
2.420e-17
As a result of variations like these are unimportant, I test that the computed outcomes are near the verified outcomes as an alternative of requiring that they be precisely the identical. (You may not see these variations for those who run this instance in your 8-processor machine. The variations rely on what else your pc is doing once you run the instance.)
Writing a certification script
I routinely use the 4 following methods to jot down a certification script:
- I test that my command reproduces outcomes that I’ve beforehand verified.
- I test that my command produces outcomes which can be near these produced by a collection of hand calculations.
- I test my command towards itself.
- I test that my command produces outcomes sufficiently shut to a different Stata command.
Certifying my command towards beforehand verified outcomes
Think about the outcomes produced by mypoisson5.ado displayed in instance 5.
Instance 5: mypoisson5 outcomes
. clear all
. use accident3
. mypoisson5 accidents cvalue children visitors
Iteration 0: f(p) = -851.18669
Iteration 1: f(p) = -556.66855
Iteration 2: f(p) = -555.81731
Iteration 3: f(p) = -555.81538
Iteration 4: f(p) = -555.81538
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .0706484 -9.28 0.000 -.7943553 -.5174188
children | -1.009017 .0807961 -12.49 0.000 -1.167374 -.8506596
visitors | .1467115 .0313762 4.68 0.000 .0852153 .2082076
_cons | .5743541 .2839515 2.02 0.043 .0178194 1.130889
------------------------------------------------------------------------------
The outcomes displayed in instance 5 are saved in e(). I beforehand verified that these outcomes are right by evaluating the higher-precision outcomes displayed in instance 6 with different outcomes.
Instance 6: mypoisson5 e() outcomes
. ereturn record
scalars:
e(N) = 505
e(rank) = 4
macros:
e(cmd) : "mypoisson5"
e(predict) : "mypoisson5_p"
e(properties) : "b V"
matrices:
e(b) : 1 x 4
e(V) : 4 x 4
features:
e(pattern)
. matrix record e(b), format(%16.15g)
e(b)[1,4]
cvalue children visitors _cons
y1 -.65588706902223 -1.0090169724739 .1467114650851 .57435412474423
I might use the outcomes from instance 6 to create a certification script like test1.do.
Code block 1: test1.do
clear all
use accident3
mypoisson5 accidents cvalue children visitors
matrix b1 = e(b)
matrix btrue = (-.65588706902223, -1.0090169724739, ///
.1467114650851, .57435412474423)
show mreldif(b1, btrue)
assert mreldif(b1, btrue) < 1e-14
Operating test1.do produces
Instance 7: test1
. do test1
. clear all
. use accident3
. mypoisson5 accidents cvalue children visitors
Iteration 0: f(p) = -851.18669
Iteration 1: f(p) = -556.66855
Iteration 2: f(p) = -555.81731
Iteration 3: f(p) = -555.81538
Iteration 4: f(p) = -555.81538
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .0706484 -9.28 0.000 -.7943553 -.5174188
children | -1.009017 .0807961 -12.49 0.000 -1.167374 -.8506596
visitors | .1467115 .0313762 4.68 0.000 .0852153 .2082076
_cons | .5743541 .2839515 2.02 0.043 .0178194 1.130889
------------------------------------------------------------------------------
. matrix b1 = e(b)
. matrix btrue = (-.65588706902223, -1.0090169724739, ///
> .1467114650851, .57435412474423)
. show mreldif(b1, btrue)
6.742e-15
. assert mreldif(b1, btrue) < 1e-14
.
.
finish of do-file
Word the method. After verifing the outcomes produced by mypoisson5, I write a certification script to make sure that mypoisson5 will all the time produce roughly these numbers. Following this course of protects me from by chance inflicting my command to provide incorrect outcomes as I make it “higher” or sooner. Don’t underestimate the significance of this safety. Placing bugs into your calculations as you try to enhance your command is remarkably simple. This course of additionally paperwork that I’ve checked this explicit case. If somebody claims to have a program that differs from mine on this case, I can ask that individual to compute the outcomes for this instance during which I do know that my command works. This request virtually all the time yields a dialogue during which that individual debugs his or her personal program in order that it produces my verified outcomes.
Right here I copied and pasted numbers from a log file right into a do-file to create test1.do. The copy-paste methodology is error–inclined, tedious, and needs to be averted. The mkassert command solves this drawback. mkassert creates assert instructions that certify outcomes saved in e(), r(), the dataset, or different Stata objects. I exploit it on a regular basis.
I start writing a certification script utilizing mkassert with code like that in test2.do, whose output seems in instance 8.
Code block 2: test2.do
clear all
use accident3
mypoisson5 accidents cvalue children visitors
mkassert eclass, mtol(1e-12) saving(test3.do, exchange)
Instance 8: Utilizing mkassert
. do test2
. clear all
. use accident3
. mypoisson5 accidents cvalue children visitors
Iteration 0: f(p) = -851.18669
Iteration 1: f(p) = -556.66855
Iteration 2: f(p) = -555.81731
Iteration 3: f(p) = -555.81538
Iteration 4: f(p) = -555.81538
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .0706484 -9.28 0.000 -.7943553 -.5174188
children | -1.009017 .0807961 -12.49 0.000 -1.167374 -.8506596
visitors | .1467115 .0313762 4.68 0.000 .0852153 .2082076
_cons | .5743541 .2839515 2.02 0.043 .0178194 1.130889
------------------------------------------------------------------------------
. mkassert eclass, mtol(1e-12) saving(test3.do, exchange)
.
finish of do-file
test2.do produces outcomes that I beforehand verified and makes use of mkassert to jot down an assert command for each consequence saved in e() within the file test3.do, which is in code block 2.
Code block 3: test3.do
assert `"`e(cmd)'"' == `"mypoisson5"'
assert `"`e(predict)'"' == `"mypoisson5_p"'
assert `"`e(properties)'"' == `"b V"'
assert e(N) == 505
assert e(rank) == 4
qui {
mat T_b = J(1,4,0)
mat T_b[1,1] = -.6558870690222316
mat T_b[1,2] = -1.009016972473914
mat T_b[1,3] = .1467114650851019
mat T_b[1,4] = .5743541247442324
}
matrix C_b = e(b)
assert mreldif( C_b , T_b ) < 1e-12
_assert_streq `"`: rowfullnames C_b'"' `"y1"'
_assert_streq `"`: colfullnames C_b'"' `"cvalue children visitors _cons"'
mat drop C_b T_b
qui {
mat T_V = J(4,4,0)
mat T_V[1,1] = .0049911902167341
mat T_V[1,2] = .0002953642487161
mat T_V[1,3] = -.0000506909358346
mat T_V[1,4] = -.0089523155601508
mat T_V[2,1] = .0002953642487161
mat T_V[2,2] = .0065280055261688
mat T_V[2,3] = .0002050149836939
mat T_V[2,4] = -.0054776138886792
mat T_V[3,1] = -.0000506909358346
mat T_V[3,2] = .0002050149836939
mat T_V[3,3] = .0009844631577381
mat T_V[3,4] = -.0075052131640854
mat T_V[4,1] = -.0089523155601508
mat T_V[4,2] = -.0054776138886792
mat T_V[4,3] = -.0075052131640854
mat T_V[4,4] = .0806284655627814
}
matrix C_V = e(V)
assert mreldif( C_V , T_V ) < 1e-12
_assert_streq `"`: rowfullnames C_V'"' `"cvalue children visitors _cons"'
_assert_streq `"`: colfullnames C_V'"' `"cvalue children visitors _cons"'
mat drop C_V T_V
Every assert command checks that what’s presently in e() is sufficiently near the corresponding worth saved in e() by mypoisson5. Strains 2–4 test the native macros, traces 6–7 test the scalars, traces 9–20 test e(b), and features 22–45 test e(V).
I create the script test4.do, which checks this case by changing the mkassert command in test2.do with the assert instructions it created in test3.do; see code block 4.
Code block 4: test4.do
// Take a look at case 1
clear all
use accident3
mypoisson5 accidents cvalue children visitors
assert `"`e(cmd)'"' == `"mypoisson5"'
assert `"`e(predict)'"' == `"mypoisson5_p"'
assert `"`e(properties)'"' == `"b V"'
assert e(N) == 505
assert e(rank) == 4
qui {
mat T_b = J(1,4,0)
mat T_b[1,1] = -.6558870690222316
mat T_b[1,2] = -1.009016972473914
mat T_b[1,3] = .1467114650851019
mat T_b[1,4] = .5743541247442324
}
matrix C_b = e(b)
assert mreldif( C_b , T_b ) < 1e-12
_assert_streq `"`: rowfullnames C_b'"' `"y1"'
_assert_streq `"`: colfullnames C_b'"' `"cvalue children visitors _cons"'
mat drop C_b T_b
qui {
mat T_V = J(4,4,0)
mat T_V[1,1] = .0049911902167341
mat T_V[1,2] = .0002953642487161
mat T_V[1,3] = -.0000506909358346
mat T_V[1,4] = -.0089523155601508
mat T_V[2,1] = .0002953642487161
mat T_V[2,2] = .0065280055261688
mat T_V[2,3] = .0002050149836939
mat T_V[2,4] = -.0054776138886792
mat T_V[3,1] = -.0000506909358346
mat T_V[3,2] = .0002050149836939
mat T_V[3,3] = .0009844631577381
mat T_V[3,4] = -.0075052131640854
mat T_V[4,1] = -.0089523155601508
mat T_V[4,2] = -.0054776138886792
mat T_V[4,3] = -.0075052131640854
mat T_V[4,4] = .0806284655627814
}
matrix C_V = e(V)
assert mreldif( C_V , T_V ) < 1e-12
_assert_streq `"`: rowfullnames C_V'"' `"cvalue children visitors _cons"'
_assert_streq `"`: colfullnames C_V'"' `"cvalue children visitors _cons"'
mat drop C_V T_V
Each time I run test4.do, it checks that mypoisson5 produces right outcomes for this one case. The extra circumstances that I confirm and certify, the extra sure I’m that my command works.
I summarize this necessary course of beneath.
- I write a do-file, right here referred to as test2.do, that produces outcomes for a case during which I’ve verified that my command produces right outcomes.
- On the finish of test2.do, I exploit mkassert to create one other do-file, right here referred to as test3.do, that accommodates assert instructions for every consequence that my command saved in e().
- I exchange the mkassert command in test2.do with the instructions it created in test3.do to create the certification script, right here referred to as test4.do.
This methodology assumes that I’ve already verified that my command produces right outcomes for a particular instance. The frequent case of verification by simulation makes this methodology extra relevant than you may suppose.
Certifying my command towards hand-calculated outcomes
I can virtually all the time discover one other strategy to compute estimation leads to Stata that needs to be numerically equal. Within the Poisson–regression case at hand, I can use gmm. As mentioned by Cameron and Trivedi (2005) and Wooldridge (2010), Poisson regression finds the (widehat{betab}) that solves the rating equations,
$$
sum_{i=1}^N left[y_i – exp(xb_iwidehat{betab})right]xb_i = {bf 0}
$$
We confirmed the way to use the gmm for related issues in Understanding the generalized methodology of moments (GMM): A easy instance. In instance 9, I exploit gmm, and I exploit assert to test that the purpose estimates produced by gmm and mypoisson5 are sufficiently shut.
Instance 9: Utilizing gmm to certify mypoisson5
. clear all
. use accident3
. mypoisson5 accidents cvalue children visitors
Iteration 0: f(p) = -851.18669
Iteration 1: f(p) = -556.66855
Iteration 2: f(p) = -555.81731
Iteration 3: f(p) = -555.81538
Iteration 4: f(p) = -555.81538
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .0706484 -9.28 0.000 -.7943553 -.5174188
children | -1.009017 .0807961 -12.49 0.000 -1.167374 -.8506596
visitors | .1467115 .0313762 4.68 0.000 .0852153 .2082076
_cons | .5743541 .2839515 2.02 0.043 .0178194 1.130889
------------------------------------------------------------------------------
. matrix b1 = e(b)
. gmm (accidents - exp({xb:cvalue children visitors _cons})), ///
> devices(cvalue children visitors) onestep
Step 1
Iteration 0: GMM criterion Q(b) = .57041592
Iteration 1: GMM criterion Q(b) = .01710408
Iteration 2: GMM criterion Q(b) = .00015313
Iteration 3: GMM criterion Q(b) = 2.190e-08
Iteration 4: GMM criterion Q(b) = 3.362e-16
word: mannequin is strictly recognized
GMM estimation
Variety of parameters = 4
Variety of moments = 4
Preliminary weight matrix: Unadjusted Variety of obs = 505
------------------------------------------------------------------------------
| Strong
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .1094934 -5.99 0.000 -.8704901 -.441284
children | -1.009017 .1884791 -5.35 0.000 -1.378429 -.6396047
visitors | .1467115 .0923401 1.59 0.112 -.0342718 .3276947
_cons | .5743542 .6039059 0.95 0.342 -.6092797 1.757988
------------------------------------------------------------------------------
Devices for equation 1: cvalue children visitors _cons
. matrix b2 = e(b)
. show mreldif(b1, b2)
5.554e-08
. assert mreldif(b1, b2) <1e-7
I used a weak tolerance when evaluating the 2 vectors of level estimates as a result of the instructions use completely different algorithms to search out their options. If I diminished the convergence tolerance in every command, the options could be nearer to one another.
For an actual certification script, I might additionally test all the pieces else saved in e() by mypoisson5 towards a price computed by gmm. I skip these particulars to current different strategies.
Certifying my command towards itself
Nearly all estimation instructions settle for if or in pattern restrictions, and these restrictions can often be examined by evaluating different outcomes produced by the identical command. Instance 10 offers an instance.
Instance 10: Testing a command towards itself
. clear all
. use accident3
. mypoisson5 accidents cvalue children visitors if cvalue <=3
Iteration 0: f(p) = -712.62548
Iteration 1: f(p) = -540.56297
Iteration 2: f(p) = -529.54572
Iteration 3: f(p) = -529.44627
Iteration 4: f(p) = -529.44618
Iteration 5: f(p) = -529.44618
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.3646368 .0872777 -4.18 0.000 -.5356979 -.1935756
children | -.9874777 .0805708 -12.26 0.000 -1.145394 -.8295618
visitors | .1488243 .0317338 4.69 0.000 .0866272 .2110214
_cons | .1081705 .3015328 0.36 0.720 -.4828229 .6991638
------------------------------------------------------------------------------
. matrix b1 = e(b)
. maintain if cvalue <=3
(121 observations deleted)
. mypoisson5 accidents cvalue children visitors
Iteration 0: f(p) = -712.62548
Iteration 1: f(p) = -540.56297
Iteration 2: f(p) = -529.54572
Iteration 3: f(p) = -529.44627
Iteration 4: f(p) = -529.44618
Iteration 5: f(p) = -529.44618
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.3646368 .0872777 -4.18 0.000 -.5356979 -.1935756
children | -.9874777 .0805708 -12.26 0.000 -1.145394 -.8295618
visitors | .1488243 .0317338 4.69 0.000 .0866272 .2110214
_cons | .1081705 .3015328 0.36 0.720 -.4828229 .6991638
------------------------------------------------------------------------------
. matrix b2 = e(b)
. show mreldif(b1, b2)
0
. assert mreldif(b1, b2) <1e-14
I start by storing the purpose estimates obtained from the pattern during which cvalue<=3 in b1. Subsequent, I maintain solely these observations within the pattern and use mypoisson5 with out an if restriction to compute the purpose estimates saved in b2. Lastly, I assert that b1 and b2 are sufficiently shut. On this case, the outcomes are precisely the identical, however I solely check that they’re shut as a result of I mustn’t depend on this equality. (I’m utilizing Stata/MP, and different jobs on my pc might change the variety of processors I successfully have, which may trigger the outcomes to vary barely.)
A similar course of works for testing in restrictions and integer weights.
Certifying my command towards one other Stata command
Generally constraining a parameter within the new estimator produces the identical outcomes as one other estimator already carried out in Stata. For instance, a random-effects estimator might cut back to a cross-sectional estimator when the variance of the random-effect is constrained to zero.
Within the case at hand, I might test that my command produces the identical values as poisson, as proven in instance 11.
Instance 11: Certifying towards an present command
. clear all
. use accident3
. mypoisson5 accidents cvalue children visitors
Iteration 0: f(p) = -851.18669
Iteration 1: f(p) = -556.66855
Iteration 2: f(p) = -555.81731
Iteration 3: f(p) = -555.81538
Iteration 4: f(p) = -555.81538
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .0706484 -9.28 0.000 -.7943553 -.5174188
children | -1.009017 .0807961 -12.49 0.000 -1.167374 -.8506596
visitors | .1467115 .0313762 4.68 0.000 .0852153 .2082076
_cons | .5743541 .2839515 2.02 0.043 .0178194 1.130889
------------------------------------------------------------------------------
. matrix b1 = e(b)
. poisson accidents cvalue children visitors
Iteration 0: log chance = -555.86605
Iteration 1: log chance = -555.8154
Iteration 2: log chance = -555.81538
Poisson regression Variety of obs = 505
LR chi2(3) = 340.20
Prob > chi2 = 0.0000
Log chance = -555.81538 Pseudo R2 = 0.2343
------------------------------------------------------------------------------
accidents | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cvalue | -.6558871 .0706484 -9.28 0.000 -.7943553 -.5174188
children | -1.009017 .0807961 -12.49 0.000 -1.167374 -.8506594
visitors | .1467115 .0313762 4.68 0.000 .0852153 .2082076
_cons | .574354 .2839515 2.02 0.043 .0178193 1.130889
------------------------------------------------------------------------------
. matrix b2 = e(b)
. show mreldif(b1, b2)
1.081e-07
. assert mreldif(b1, b2) <1e-6
Finished and undone
I offered some methods that I exploit to jot down certification scripts. An actual certification script would cowl many extra circumstances. Within the subsequent publish, I focus on utilizing and creating Mata libraries.
References
Cameron, A. C., and P. Okay. Trivedi. 2005. Microeconometrics: Strategies and Functions. Cambridge: Cambridge College Press.
Gould, W. 2001. Statistical software program certification. Stata Journal 1: 29–50.
Wooldridge, J. M. 2010. Econometric Evaluation of Cross Part and Panel Knowledge. 2nd ed. Cambridge, Massachusetts: MIT Press.