¿Regresión para un modelo de forma ?

Tengo un conjunto de datos que son estadísticas de un foro de discusión web. Estoy viendo la distribución de la cantidad de respuestas que se espera que tenga un tema. En particular, he creado un conjunto de datos que tiene una lista de conteos de respuestas de temas, y luego el conteo de temas que tienen ese número de respuestas.

"num_replies","count"
0,627568
1,156371
2,151670
3,79094
4,59473
5,39895
6,30947
7,23329
8,18726

Si trazo el conjunto de datos en un diagrama de log-log, obtengo lo que es básicamente una línea recta:

(Esta es una distribución Zipfian ). Wikipedia me dice que las líneas rectas en las gráficas log-log implican una función que puede ser modelada por un monomio de la forma . Y, de hecho, he analizado dicha función:$y = ax^k$

lines(data$num_replies, 480000 * data$num_replies ^ -1.62, col="green")

Mis globos oculares obviamente no son tan precisos como R. Entonces, ¿cómo puedo hacer que R se ajuste a los parámetros de este modelo para mí con mayor precisión? Intenté la regresión polinómica, pero no creo que R intente ajustar el exponente como parámetro: ¿cuál es el nombre apropiado para el modelo que quiero?

Editar: Gracias por las respuestas a todos. Como se sugirió, ahora he ajustado un modelo lineal contra los registros de los datos de entrada, usando esta receta:

data <- read.csv(file="result.txt")

# Avoid taking the log of zero:
data$num_replies = data$num_replies + 1

plot(data$num_replies, data$count, log="xy", cex=0.8)

# Fit just the first 100 points in the series:
model <- lm(log(data$count[1:100]) ~ log(data$num_replies[1:100]))

points(data$num_replies, round(exp(coef(model)[1] + coef(model)[2] * log(data$num_replies))), 
       col="red")

El resultado es este, mostrando el modelo en rojo:

Eso parece una buena aproximación para mis propósitos.

Si luego uso este modelo Zipfian (alfa = 1.703164) junto con un generador de números aleatorios para generar el mismo número total de temas (1400930) que contenía el conjunto de datos medido original (usando este código C que encontré en la web ), el resultado se ve me gusta:

Los puntos medidos están en negro, los generados aleatoriamente según el modelo están en rojo.

Creo que esto muestra que la varianza simple creada al generar aleatoriamente estos 1400930 puntos es una buena explicación para la forma del gráfico original.

Si está interesado en jugar con los datos sin procesar usted mismo, los he publicado aquí .

Su ejemplo es muy bueno porque señala claramente problemas recurrentes con dichos datos.

Dos nombres comunes son función de poder y ley de poder. En biología, y en algunos otros campos, la gente suele hablar de alometría, especialmente cuando relaciona medidas de tamaño. En física y en otros campos, la gente habla de leyes de escala.

No consideraría monomial como un buen término aquí, ya que lo asocio con poderes enteros. Por la misma razón, esto no se considera mejor como un caso especial de un polinomio.

Los problemas de ajustar una ley de poder a la cola de una distribución se transforman en problemas de ajustar una ley de poder a la relación entre dos variables diferentes.

La forma más fácil de ajustar una ley de potencia es tomar logaritmos de ambas variables y luego ajustar una línea recta mediante regresión. Hay muchas objeciones a esto siempre que ambas variables estén sujetas a error, como es común. El ejemplo aquí es un caso puntual ya que ambas variables (y ninguna) podrían considerarse como respuesta (variable dependiente). Ese argumento conduce a un método de ajuste más simétrico.

Además, siempre está la cuestión de los supuestos sobre la estructura del error. Una vez más, el ejemplo aquí es un ejemplo: los errores son claramente heteroscedasticos. Eso sugiere algo más como mínimos cuadrados ponderados.

Una excelente revisión es http://www.ncbi.nlm.nih.gov/pubmed/16573844

Otro problema es que las personas a menudo identifican leyes de poder solo en un rango de sus datos. Luego, las preguntas se vuelven científicas y estadísticas, hasta determinar si las leyes de poder son solo una ilusión o un pasatiempo aficionado de moda. Gran parte de la discusión surge bajo los encabezados de comportamiento fractal y sin escala, con una discusión asociada que va desde la física hasta la metafísica. En su ejemplo específico, una pequeña curvatura parece evidente.

Los entusiastas de las leyes de poder no siempre son igualados por los escépticos, porque los entusiastas publican más que los escépticos. Sugeriría que un diagrama de dispersión en escalas logarítmicas, aunque es un diagrama natural y excelente que es esencial, debe ir acompañado de gráficos residuales de algún tipo para verificar las desviaciones de la forma de la función de potencia.

Si supone que una potencia es un buen modelo para ajustarse, puede usarla log(y) ~ log(x)como modelo y ajustar una regresión lineal usando lm():

Prueba esto:

# Generate some data
set.seed(42)

x <- seq(1, 10, 1)

a = 10
b = 2
scatt <- rnorm(10, sd = 0.2)


dat <- data.frame(
  x = x,
  y = a*x^(-b) + scatt
)

Ajustar un modelo:

# Fit a model
model <- lm(log(y) ~ log(x) + 1, data = dat) 
summary(model)

pred <- data.frame(
  x = dat$x,
  p = exp(predict(model, dat))
)

Ahora crea una trama:

# Create a plot
library(ggplot2)
ggplot() +
  geom_point(data = dat, aes(x=x, y=y)) +
  geom_line(data = pred, aes(x=x, y=p), col = "red")