Fitting of Negative Binomial

cat("Fitting of Negative Bnomial");
cat("Let X~NB(x;k,p)");

x = c(0,1,2,3,4,5);
f = c(450,132,41,22,3,2);
fx = fx;
fx2= f(x^2);

data = data.frame(x,f,fx,fx2);data

mean = sum(fx)/sum(f);
var = (sum(fx2)/sum(f))-(mean)^2;
N = sum(f);
sumfx2 = sum(fx2);
mean2 = mean*mean;
cat(paste("N = ", N));
cat(paste("mean = ", mean));
cat(paste("var = ", var));
cat(paste("sumfx2 = ", sumfx2));
cat(paste("mean^2 = ", mean2));

p = mean/var;
q = 1 - p;
k = mean*p/q;

cat(paste("p=",p,"q=",q,"k=",k));

cat("Reccurence rel");
cat("a = x, b = f, c = ( (x + k)/(x+1) )q, d = P[X=x+1], e = Expected freq = Np(X=x)");

c = ((x+k)/(x+1))*q;
d = dnbinom(x,k,p);
ef = sum(f)*d;

data = data.frame(x,f,c,d,ef);data

psum = sum(d);
efsum = sum(ef);
cat(paste("psum = ", psum));
cat(paste("efsum = ", efsum));