Thanks for the feedback. I'll check with Bobby.
For those who don't know my purpose in asking, I'm interested in programming
a PIG (or Holla) model in Stata, which is a Poisson-Inverse Gaussian mixture
model. I did it using MathStatica 1.5, which is a statistical add-on package
to Mathematica 5.2. But since a type 2 Bessel function is involved in both
the PMF and LL functions, it may be difficult to program it using Stata
facilities. Hence my query about any good approximations - or even a programmed
Bessel type 2 function itself. Since the Bessel is based on the Gamma, I thought
perhaps someone might have found a shortcut to the Bessel using Stata's
gamma - or a related - function.
I found an approximation to the Bessel Type 2 in an old journal article.
Since it is an approximation anyhow, I further simplified it by calculating pi/2
and pi/4 to their numeric values (to 4 decimal points), which are 1.5708 and
0.7854 respectively. The actual approximation I found was given as:
J_n(x) ~ sqrt(2/(pi*x)) * sin(x-(2*n+1)*(pi/4).
Kit Baum provided me with a reference to Fortran source code for the
function. I'll take a look at that as well. It just may be that there will be too
much work involved to program the PIG at this time. However, it is has useful
properties for modeling the response as the values of the counts slowly
decrease with their increasing value.
Joe Hilbe
=========================================================
Date: Sat, 1 Oct 2005 17:43:21 +0100
From: "Nick Cox" <[email protected]>
Subject: st: RE: RE: Bessel
Bobby Gutierrez needed some Bessel function
for some purpose and wrote an undocumented helper
program. I can't remember which and can't find
the code.
I once wrote a rather specific Stata program for I_0,
not what you specify here.
1.5708 means _pi/2 here.
Odds are you need to transliterate formulae
from Abramowitz and Stegun's handbook to Stata or Mata.
Nick
[email protected]
[email protected]
> To be more specific, has anyone made a Bessel of the 2nd
> kind function in
> Stata? An approximation is sqrt(1.5708/x) *
> sin(x-0.7854(2y+1)) for values
> closer to the max values and
> -[ (2^n)(n-1)! /pi ] * x^(-n) for values closer to 0. If
> someone has a
> better approximation, or the "exact" function in Stata, I'd
> love to see it.
*
* 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/
