Update:

@J.M.'stechnicaldifficulties noted in a comment that the distribution could be written as follows:

BETDistribution[x_, c_] := TransformedDistribution[r1 (1 + r2), 
  {r1 \[Distributed] BernoulliDistribution[1 - (1 - x)/(1 + (c - 1) x)],
   r2 \[Distributed] GeometricDistribution[1 - x]}, 
   Assumptions -> c >= 1 && 0 < x < 1]

Then this allows RandomVariate to work properly:

SeedRandom[12345];
data = RandomVariate[BETDistribution[1/2, 5], 1000];

So no need for writing one's own functions to obtain random samples.

But there is one unforeseen downside: FindDistributionParameters is much, much slower with this definition of BETDistribution. With the above data and the newer definition of BETDistribution we have the following:

AbsoluteTiming[FindDistributionParameters[data, BETDistribution[x, c]]]
(* {22.7427, {x -> 0.505552, c -> 5.37284}} *)

With the other definition we have

BETDistribution[x_, c_] := ProbabilityDistribution[Piecewise[{{(1 - x)/(1 + (c - 1) x),
  k == 0}}, c (1 - x) (x^k)/(1 + (c - 1) x)], {k, 0, ∞, 1}, 
  Assumptions -> c >= 1 && 0 < x < 1]

AbsoluteTiming[FindDistributionParameters[data, BETDistribution[x, c]]]
(* {0.0748486, {c -> 5.37284, x -> 0.505552}} *)

That's 300 times longer with the TransformedDistribution. (The Rolling Stones said it long ago: "You can't always get what you want.")

@J.M.'stechnicaldifficulties answer showed how to use Piecewise to obtain the desired definition which then allows FindDistributionParameters to work. But the question of generating random samples from this distribution still remains.

In Mathematica 12.1

RandomVariate[BETDistribution[1/2, 5], 10]

returns

Fortunately in this case it is relatively easy and quick to generate a large random sample. We separate the random selection of 0's and non-0's. First a Bernoulli random number is selected with probability $1 - Pr[0] = 1 - (1 - x)/(1 + (-1 + c) x)$. If that random number is zero, then 0 is selected. If not, then it turns out that the random variable $Z|Z>0$ has the same distribution of 1 plus a Geometric random variable with parameter 1 - x. Such a function can be written as

rvBET[x_, c_, nSamples_] := Module[{z1, z2},
  z1 = RandomVariate[BernoulliDistribution[1 - (1 - x)/(1 + (c - 1) x)], nSamples];
  z2 = 1 + RandomVariate[GeometricDistribution[1 - x], nSamples];
  z1*z2
  ]

As a partial check on this consider generating a large amount of data with known parameters and then attempt to estimate the parameters:

SeedRandom[12345];
data = rvBET[1/4, 5, 100000];
FindDistributionParameters[data, BETDistribution[x, c]]
(* {c -> 4.9875, x -> 0.251256} *)

@J.M.'stechnicaldifficulties answer showed how to use Piecewise to obtain the desired definition which then allows FindDistributionParameters to work. But the question of generating random samples from this distribution still remains.

In Mathematica 12.1

RandomVariate[BETDistribution[1/2, 5], 10]

returns

Fortunately in this case it is relatively easy and quick to generate a large random sample. We separate the random selection of 0's and non-0's. First a Bernoulli random number is selected with probability $1 - Pr[0] = 1 - (1 - x)/(1 + (-1 + c) x)$. If that random number is zero, then 0 is selected. If not, then it turns out that the random variable $Z|Z>0$ (where $Z\sim \text{BETDistribution}(x,c)$) has the same distribution of 1 plus a Geometric random variable with parameter 1 - x. Such a function can be written as

rvBET[x_, c_, nSamples_] := Module[{z1, z2},
  z1 = RandomVariate[BernoulliDistribution[1 - (1 - x)/(1 + (c - 1) x)], nSamples];
  z2 = 1 + RandomVariate[GeometricDistribution[1 - x], nSamples];
  z1*z2
  ]

As a partial check on this consider generating a large amount of data with known parameters and then attempt to estimate the parameters:

SeedRandom[12345];
data = rvBET[1/4, 5, 100000];
FindDistributionParameters[data, BETDistribution[x, c]]
(* {c -> 4.9875, x -> 0.251256} *)