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/
