Here is the H-G probability density,  written as a function of µ, which is cos θ.

pdf[µ] = 1/2 (1 - g^2) / (1 + g^2 - 2 g µ)^(3/2)

Now we can get the cumulative distribution by integrating with respect to μ.

cdf[µ'] = Integrate[pdf[µ], {µ, -1, µ'}, Assumptions->{-1 < g < 1, -1 < µ' < 1, g µ' ≠ 0}]

The sampling procedure is to solve for cdf equal to ξ where ξ is a uniform random number between 0 and 1.

sample[ξ] = µ' /. First[Solve[cdf[µ'] == ξ, µ']]

This is the sampling function as I have it in my notes:

mysample[ξ] = (-((1-g^2)/(1 - g + 2 g ξ))^2 + (1 + g^2))/(2g)

Verify that it is the same as the one just derived.

Together[mysample[ξ] == sample[ξ]]

Here is a simplified version that may be better for evaluation:

Simplify[sample[ξ]]