More general than the ! gamma J function, we can use:
gma=: 4 : 0
'a b'=.x
(((y%b)^<:a)*^-y%b)%b*!<:a
)
where a is the shape parameter and b is the scale parameter.
I am looking for a Gamma distribution randon number generator in J.
Wikipedia suggests:
We provide an algorithm without proof. It is an instance of the
acceptance-rejection
method <http://en.wikipedia.org/wiki/Rejection_sampling>:
1. Let *m* be 1.
2. Generate *V*3*m* − 2, *V*3*m* − 1 and *V*3*m* ― independent uniformly
distributed on (0, 1] variables.
3. If [image: V_{3m - 2} \le v_0], where [image: v_0 = \frac e {e +
\delta}], then go to step 4, else go to step 5.
4. Let [image: \xi_m = V_{3m - 1}^{1 / \delta}, \ \eta_m = V_{3m} \xi _m^
{\delta - 1}]. Go to step 6.
5. Let [image: \xi_m = 1 - \ln {V_{3m - 1}}, \ \eta_m = V_{3m}
e^{-\xi_m}].
6. If [image: \eta_m > \xi_m^{\delta - 1} e^{-\xi_m}], then
increment *m*and go to step 2.
7. Assume ξ = ξ*m* to be the realization of Γ(δ,1)
Now, to summarize,
[image: \theta \left( \xi - \sum _{i=1} ^{[k]} {\ln U_i} \right) \sim \Gamma
(k, \theta),]
where [*k*] is the integral part of *k*, and *ξ* has been generated using
the algorithm above with δ = {*k*} (the fractional part of *k*), *Uk* and *V
l* are distributed as explained above and are all independent.
In J I tried:
getG=: 3 : 0
e=.0
a=. (-<.)y
if.(0~:a) do.
whilst.(n>(em*^<:a)*^-e) do.
v=.(/:)~?3$0
em=.e
if.(({.v)<%>:a%e.1)do.
e=.(1{.v)^%a
n=.({:v)*e^<:a
else.
e=.-.^.1{.v
n=.({:v)*^-e
end.
end.
end.
+/e,-^. ?(<.y)#0
)
This thing would be sluggish if it worked, but alas it does not. When y is
not an integer, it seems to go into an infinite loop.
If this attempt fails, i may have to revert to a more crude routine like the
following for n random numbers
smp=.0.001*(+/\(2.2,1)gma"1 0( 0.001*i.10000))I.?n#1000
Plotting the frequency distribution of these numbers should produce
something like a Gamma if n is large enough with
plot {:"1 /:~({.,#)/.~0.1*<.0.5+10*smp
In the past I enjoyed using J for a simple MonteCarlo type simulation. In
an attempt to move on to something more complex I am comparing J to
SSJ(Statistical Simulation in Java by Pierre l'Ecuyer). On the first sample
in this SSJ pachage a J version was easy and faster, but the second required
Gamma random numbers.
Any suggestions?
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
