Regresión para un modelo de forma $y=ax^k$ ?

Question

Tengo un conjunto de datos que son estadísticas de un foro de discusión de la web. Estoy buscando la distribución del número de respuestas que se espera que tenga un tema. En particular, he creado un conjunto de datos que tiene una lista de recuentos de respuestas de temas, y luego el recuento 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 gráfico logarítmico, obtengo lo que es básicamente una línea recta:

(Este es un Distribución zipfiana ). Wikipedia me dice que las líneas rectas en los gráficos logarítmicos implican una función que puede ser modelada por un monomio de la forma $y = ax^k$ . Y de hecho he ojeado tal función:

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

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

Edición: Gracias a todos por las respuestas. Como se sugirió, ahora he ajustado un modelo lineal contra los registros de los datos de entrada, utilizando 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:

Parece una buena aproximación para mis propósitos.

Si luego utilizo este modelo Zipfiano (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 originalmente (utilizando este código C que encontré en la web ), el resultado es el siguiente:

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

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

Si estás interesado en jugar con los datos en bruto por ti mismo, tengo lo publicó aquí .

Answer 1

Su ejemplo es muy bueno porque señala claramente los problemas recurrentes de estos datos.

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

Yo no consideraría que el monomio es un buen término aquí, ya que lo asocio con potencias enteras. Por la misma razón no es mejor considerarlo como un caso especial de un polinomio.

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

La forma más fácil de ajustar una ley de potencia es tomar los logaritmos de ambas variables y luego ajustar una línea recta utilizando la regresión. Hay muchas objeciones a esto siempre que ambas variables estén sujetas a error, como es habitual. El ejemplo que nos ocupa es un buen ejemplo, ya que ambas variables (y ninguna) podrían considerarse como respuesta (variable dependiente). Este argumento conduce a un método de ajuste más simétrico.

Además, siempre está la cuestión de las suposiciones sobre la estructura del error. Una vez más, el ejemplo aquí es un buen ejemplo, ya que los errores son claramente heteroscedásticos. Eso sugiere algo más parecido a los mínimos cuadrados ponderados.

Una excelente reseña es http://www.ncbi.nlm.nih.gov/pubmed/16573844

Otro problema es que la gente suele identificar las leyes de potencia sólo en algún rango de sus datos. Las cuestiones se convierten entonces en científicas y estadísticas, hasta llegar a la cuestión de si la identificación de las leyes de potencia es sólo una ilusión o un pasatiempo de aficionados de moda. Gran parte de la discusión surge bajo los epígrafes de comportamiento fractal y sin escala, con un debate asociado que va desde la física hasta la metafísica. En tu ejemplo concreto, parece evidente una pequeña curvatura.

Los entusiastas de las leyes de la energía no siempre se encuentran con los escépticos, porque los entusiastas publican más que los escépticos. Yo sugeriría que un gráfico de dispersión en escalas logarítmicas, aunque es un gráfico natural y excelente que es esencial, debería ir acompañado de gráficos residuales de algún tipo para comprobar las desviaciones de la forma de la función de potencia.

Answer 2

Si se asume que una potencia es un buen modelo para ajustar, entonces se puede utilizar log(y) ~ log(x) como su modelo, y ajuste una regresión lineal utilizando 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
)

Encaja 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 parcela:

# 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")