こんな感じになるのかな?
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