And a slightly less rough one:
program poissontt
version 9
syntax varlist(min=2 max=2), Generate(str)
confirm new var `generate'
tokenize `varlist'
args observed expected
tempvar maxprob minprob
quietly {
gen double `maxprob' = exp(`expected')
su `observed', meanonly
local max = r(max)
local max0 = r(max) - 1
forval I = 0/`max0' {
replace `maxprob' =
`maxprob' * (`expected'^`I') / exp(lnfactorial(`I')) if `observed' <= `I'
}
replace `maxprob' = 1 if `observed' == 0
replace `maxprob' = 1e-6 if `expected' < 1e-6
gen `minprob' = 0
forval I = 0/`max' {
replace `minprob' = ///
`minprob' + ((exp(-`expected')) * (`expected'^`I')) / exp(lnfactorial(`I')) ///
if `observed' <= `I'
}
gen `generate' = min(`maxprob', `minprob')
replace `generate' = 0 if `observed' == 0 & `expected' == 0
replace `generate' = max(2 * `generate', 1)
}
end
Nick
[email protected]
> Here is a very rough sketch:
>
> program poissontt
> version 9
> syntax varlist(min=2 max=2), Generate(str)
> confirm new var `generate'
> tempvar maxprob minprob
>
> quietly {
> gen double `maxprob' = exp(`expected')
>
> su `observed', meanonly
> local max = r(max)
> local max0 = r(max) - 1
>
> forval I = 0/`max0' {
> replace `maxprob' =
> `maxprob' * (`expected'^`I') / exp(lnfactorial(`I')) if
> `observed' <= `I'
> }
> replace `maxprob' = 1 if `observed' == 0
> replace `maxprob' = 1e-6 if `expected' < 1e-6
>
> gen `minprob' = 0
> forval I = 0/`max' {
> replace `minprob' = ///
> `minprob' + ((exp(-`expected')) * (`expected'^`I')) /
> exp(lnfactorial(`I')) ///
> if `observed' <= `I'
> }
>
> gen `generate' = min(`maxprob', `minprob')
> replace `generate' = 0 if `observed' == 0 &
> `expected' == 0
> replace `generate' = max(2 * `generate', 1)
> }
> end
>
> Nick
> [email protected]
>
> > A broader comment is that your approach of making this
> > byable and then looping over the data does not look
> > the best strategy to me, as calculations look as if they can
> > and should be recast as calculations on variables.
> >
>
> *
> * For searches and help try:
> * http://www.stata.com/support/faqs/res/findit.html
> * http://www.stata.com/support/statalist/faq
> * http://www.ats.ucla.edu/stat/stata/
>
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
