How to calculate Gamma functions values for Quaternions using Mathematica?

If the function to be applied to a quaternion is holomorphic/meromorphic, one way to proceed is to switch to the complex matrix representation of a quaternion (see e.g. this), and then apply MatrixFunction[] to it.

For instance,

<< Quaternions`

BlockRandom[SeedRandom[42]; (* for reproducibility *)
            qq = Quaternion @@ RandomVariate[NormalDistribution[], 4]]
   Quaternion[-1.01205, 0.826338, -1.39379, 0.41692]

Tan[qq] (* elementary example *)
   Quaternion[-0.0652529, 0.50833, -0.857402, 0.256472]

Quaternion @@ Flatten[ReIm[First[MatrixFunction[Tan,
                                        {{qq[[1]] + I qq[[2]], qq[[3]] + I qq[[4]]},
                                         {-qq[[3]] + I qq[[4]], qq[[1]] - I qq[[2]]}}]]]]
   Quaternion[-0.0652529, 0.50833, -0.857402, 0.256472]

Thus, to evaluate a gamma function for quaternion argument:

Quaternion @@ Flatten[ReIm[First[MatrixFunction[Gamma,
                                        {{qq[[1]] + I qq[[2]], qq[[3]] + I qq[[4]]},
                                         {-qq[[3]] + I qq[[4]], qq[[1]] - I qq[[2]]}}]]]]
   Quaternion[-0.0568728, 0.0210734, -0.0355446, 0.0106324]

---

Now, for a warning. The reason I emphasized the requirement of analyticity is that if the function has branch cuts, then the complex matrix method might give results that are the conjugate of the expected result, due to possible differences in branch cut choice. For instance,

Log[qq] (* a consistent example *)
   Quaternion[0.670588, 1.04449, -1.76174, 0.526985]

Quaternion @@ Flatten[ReIm[First[MatrixFunction[Log,
                                        {{qq[[1]] + I qq[[2]], qq[[3]] + I qq[[4]]},
                                         {-qq[[3]] + I qq[[4]], qq[[1]] - I qq[[2]]}}]]]]
   Quaternion[0.670588, 1.04449, -1.76174, 0.526985]

q2 = Reverse[qq] (* inconsistent example *)

Log[q2]
   Quaternion[0.670588, 0.989246, -0.586497, 0.718306]

(* note that this returns the conjugate *)
Quaternion @@ Flatten[ReIm[First[MatrixFunction[Log,
                                        {{q2[[1]] + I q2[[2]], q2[[3]] + I q2[[4]]},
                                         {-q2[[3]] + I q2[[4]], q2[[1]] - I q2[[2]]}}]]]]
   Quaternion[0.670588, -0.989246, 0.586497, -0.718306]

I have decided to write another answer, because the following is more specific to the case of the gamma function. In particular, I will demonstrate how to modify the OP's original approach.

The reason the OP's method failed is that the elementary functions in Mathematica are left unevaluated for symbolic quaternion arguments. For instance,

<<Quaternions`

Exp[Quaternion[a, b, c, d]]
   E^Quaternion[a, b, c, d]

To make the OP's proposal work, one will need to manually split the quaternion into its scalar and vector parts. Using e.g. the formula here, we can write the following:

ig[Quaternion[a_, b_, c_, d_], x_] := 
   With[{l = Log[x], vs = Sqrt[{b, c, d}.{b, c, d}]},
        (x^(a - 1) Exp[-x]) Prepend[l {b, c, d} Sinc[l vs], Cos[l vs]]]

Then, we can do this:

qq = Quaternion[1, 2, 3, 0];

Quaternion @@ Integrate[ig[qq, x], {x, 0, Infinity}]
   Quaternion[1/2 (Gamma[1 - I Sqrt[13]] + Gamma[1 + I Sqrt[13]]),
              (I (Gamma[1 - I Sqrt[13]] - Gamma[1 + I Sqrt[13]]))/Sqrt[13],
              (3 I (Gamma[1 - I Sqrt[13]] - Gamma[1 + I Sqrt[13]]))/(2 Sqrt[13]), 0]

N[%] // Chop
   Quaternion[-0.003441450173902947, 0.008960473460683197, 0.013440710191024793, 0]

Compare with the matrix function approach:

Quaternion @@ Flatten[ReIm[First[MatrixFunction[Gamma,
                        {{qq[[1]] + I qq[[2]], qq[[3]] + I qq[[4]]},
                         {-qq[[3]] + I qq[[4]], qq[[1]] - I qq[[2]]}}]]]] // FullSimplify
   Quaternion[1/2 (Gamma[1 - I Sqrt[13]] + Gamma[1 + I Sqrt[13]]),
              (I (Gamma[1 - I Sqrt[13]] - Gamma[1 + I Sqrt[13]]))/Sqrt[13],
              (3 I (Gamma[1 - I Sqrt[13]] - Gamma[1 + I Sqrt[13]]))/(2 Sqrt[13]), 0]

and we see the two methods give consistent results.

With my GeometricAlgebra paclet you can work with any geometric number and apply any MatrixFunction supported function on them with MultivectorFunction , that effectively just converts numbers to matrices and back. Quaternions can be simply expressed as GeometricAlgebra[0,2]:

PacletInstall["https://wolfr.am/OkONsyY2"]

<< GeometricAlgebra`

MultivectorFunction[Gamma, 
  Multivector[{a, b, c, d}, {0, 2}]] // Simplify

Convert quaternion to matrix:

MultivectorMatrix[Multivector[{a, b, c, d}, {0, 2}]]["Numeric"]