In case anyone is still interested, I have managed to implement Aristizabal's formulae in Java. This is more proof-of-concept than the requested "robust" code, but it is a starting point.
/**
* Computes the point estimate of the shift offset (gamma) from the given sample. The sample array will be sorted by this method.<p>
* Cf. Aristizabal section 2.2 ff.
* @param sample {@code double[]}, will be sorted
* @return gamma point estimate
*/
public static double pointEstimateOfGammaFromSample(double[] sample) {
Arrays.sort(sample);
DoubleUnaryOperator func = x->calculatePivotalOfSortedSample(sample, x)-1.0;
double upperLimit = sample[0];
double lowerLimit = 0;
double gamma = bisect(func, lowerLimit, upperLimit);
return gamma;
}
/**
* Cf. Aristizabal's equation (2.3.1)
* @param sample {@code double[]}, should be sorted in ascending order
* @param gamma shift offset
* @return pivotal value of sample
*/
private static double calculatePivotalOfSortedSample(final double[] sample, double gamma) {
final int n=sample.length;
final int n3=n/3;
final double mid = avg(sample, gamma, n3+1, n-n3);
final double low = avg(sample, gamma, 1, n3);
final double upp = avg(sample, gamma, n-n3+1, n);
final double result = (mid-low)/(upp-mid);
return result;
}
/**
* Computes average of sample values from {@code sample[l-1]} to {@code sample[u-1]}.
* @param sample {@code double[]}, should be sorted in ascending order
* @param gamma shift offset
* @param l lower limit
* @param u upper limit
* @return average
*/
private static double avg(double[] sample, double gamma, int l, int u) {
double sum = 0.0;
for (int i=l-1;i<u;sum+=Math.log(sample[i++]-gamma));
final int n = u-l+1;
return sum/n;
}
/**
* Naive bisection implementation. Should always complete if the given values actually straddles the root.
* Will call {@link #secant(DoubleUnaryOperator, double, double)} if they do not, in which case the
* call may not complete.
* @param func Function solve for root value
* @param lowerLimit Some value for which the given function evaluates < 0
* @param upperLimit Some value for which the given function evaluates > 0
* @return x value, somewhere between the lower and upper limits, which evaluates close enough to zero
*/
private static double bisect(DoubleUnaryOperator func, double lowerLimit, double upperLimit) {
final double eps = 0.000001;
double low=lowerLimit;
double valAtLow = func.applyAsDouble(low);
double upp=upperLimit;
double valAtUpp = func.applyAsDouble(upp);
if (valAtLow*valAtLow>0) {
// Switch to secant method
return secant(func, lowerLimit, upperLimit);
}
System.out.printf("bisect %f@%f -- %f@%f%n", valAtLow, low, valAtUpp, upp);
double mid;
while(true) {
mid = (upp+low)/2;
if (Math.abs(upp-low)/low<eps)
break;
double val = func.applyAsDouble(mid);
if (Math.abs(val)<eps)
break;
if (val<0)
low=mid;
else
upp=mid;
}
return mid;
}
/**
* Naive secant root solver implementation. May not complete if root not found.
* @param f Function solve for root value
* @param a Some value for which the given function evaluates
* @param b Some value for which the given function evaluates
* @return x value which evaluates close enough to zero
*/
static double secant(final DoubleUnaryOperator f, double a, double b) {
double fa = f.applyAsDouble(a);
if (fa==0)
return a;
double fb = f.applyAsDouble(b);
if (fb==0)
return b;
System.out.printf("secant %f@%f -- %f@%f%n", fa, a, fb, b);
if (fa*fb<0) {
return bisect(f, a, b);
}
while ( abs(b-a) > abs(0.00001*a) ) {
final double m = (a+b)/2;
final double k = (fb-fa)/(b-a);
final double fm = f.applyAsDouble(m);
final double x = m-fm/k;
if (Math.abs(fa)<Math.abs(fb)) {
// f(a)<f(b); Choose x and a
b=x;
fb=f.applyAsDouble(b);
} else {
// f(a)>=f(b); Choose x and b
a=x;
fa=f.applyAsDouble(a);
}
if (fa==0)
return a;
if (fb==0)
return b;
if (fa*fb<0) {
// Straddling root; switch to bisect method
return bisect(f, a, b);
}
}
return (a+b)/2;
}
The difference is what we are interested in. Both distributions are built from independent Bernoulli trials with fixed probability of success, p.
With the Binomial distribution, the random variable X is the number of successes observed in n trials. Because there are a fixed number of trials, the possible values of X are 0, 1, ..., n.
With the Negative Binomial distribution, the random variable Y is the number of trials until observed the r th success is observed. In this case, we keep increasing the number of trials until we reach r successes. The possible values of Y are r, r+1, r+2, ... with no upper bound. The Negative Binomial can also be defined in terms of the number of failures until the r th success, instead of the number of trials until the r th success. Wikipedia defines the Negative Binomial distribution in this manner.
So to summarize:
Binomial:
- Fixed number of trials (n)
- Fixed probability of success (p)
- Random variable is X = Number of successes.
- Possible values are 0 ≤ X ≤ n
Negative Binomial:
- Fixed number of successes (r)
- Fixed probability of success (p)
- Random variable is Y = Number of trials until the r th success.
- Possible values are r ≤ Y
Thanks to Ben Bolker for reminding me to mention the support of the two distributions. He answered a related question here.
Best Answer
I guess what you are looking for is the probability generating function. A derivation of the probability generating function of the binomial distribution can be found under
http://economictheoryblog.com/2012/10/21/binomial-distribution/
However, having a look at Wikipedia is nowadays always a good idea, although I have to say that the specification of the binomial could be improved.
https://en.wikipedia.org/wiki/Binomial_distribution#Specification