ダメ出し：よく分からないけど，ダメ その2

「バク＝スネッペンのゲーム」

こんな感じになるのかな?

n <- 100
limit <- 2/3
x <- runif(n+2, 0, 1)
repeat {
    xmin <- which.min(x[2:(n+1)])
    if (x[xmin+1] >= limit) break
    x[xmin:(xmin +2)] <- runif(3)
}
z <- numeric(n)
for (i in 1:n) {
    y <- 0
    xmin <- which.min(x[2:(n+1)])
    x[xmin:(xmin +2)] <- runif(3)
    repeat {
        xmin <- which.min(x[2:(n+1)])
        if (x[xmin+1] >= limit) break
        x[xmin:(xmin +2)] <- runif(3)
        y <- y+1
    }
    z[i] <- y
}
plot(sort(log(z)))

fit は，データに 0 があるとまずいのではないでしょうかね...

> z <- z[z > 0]
> summary(power.law.fit(z))
Maximum likelihood estimation

Call:
mle(minuslogl = mlogl, start = list(alpha = start))

Coefficients:
       Estimate Std. Error
alpha 0.4585079 0.05960587