def METHOD_NAME(seed, ctx, func_name, k, theta, shape):
    with nn.context_scope(ctx):
        o = F.rand_gamma(k, theta, shape, seed=seed)
    assert o.shape == tuple(shape)
    assert o.parent.name == func_name
    o.forward()
    if o.size > 10000:
        est_mu = o.d.mean()
        est_sigma = o.d.std()
    else:
        data = []
        for i in range(1000000//o.size):
            o.forward()
            data += [o.d.copy()]
        est_mu = np.mean(np.array(data))
        est_sigma = np.std(np.array(data))
    mu = k * theta  # theoretical mean
    var = k * theta * theta
    sigma = np.sqrt(var)  # theoretical std
    assert np.isclose(est_mu, mu, atol=5e-2)
    assert np.isclose(est_sigma, sigma, atol=5e-2)
    # Checking recomputation
    func_args = [k, theta, shape, seed]
    recomputation_test(rng=None, func=F.rand_gamma, vinputs=[],
                       func_args=func_args, func_kwargs={}, ctx=ctx)