Estimating the exponent from empirical data in R

   pwrdist <- function(u,...) {

      # u is vector of event counts, e.g. how many

      # crimes was a given perpetrator charged for by the police

      fx <- table(u)

      i <- as.numeric(names(fx))

      y <- rep(0,max(i))

      y[i] <- fx

      m0 <- glm(y~log(1:max(i)),family=quasipoisson())

      print(summary(m0))

      sub <- paste("s=",round(m0$coef[2],2),"lambda=",sum(u),"/",length(u))

      plot(i,fx,log="xy",xlab="x",sub=sub,ylab="counts",...)

      grid()

      lines(1:max(i),(fitted(m0)),type="l")

      return(m0)

  }