I'm having a problem and I can't figure out why the code behaves as it does. I'm simulating a clustered process in R using the package spatstat. The cluster process is the Thomas process $\mathcal{T}(\kappa,\mu,\sigma),$ thus I can simulate this process using the function rThomas.

Then, I want to estimate the parameters above using Ripley's $K$-function. This is done by using the function thomas.extK. If I fix the parameters $\kappa = 10, \ \sigma=0.1, \ \mu=6$ and generate the Thomas process $100$ times, and each time estimate the parameters, store them in a vector, then take the mean of each parametervector, I expect to get close to the fixed values. However my result is way off. I don't see any error in my code and I can't explain the discrepancy.

This is the same as if I generate $1000$ number from a normal distribution with mean $50$ and when I take the mean of all the generated points I should get close to $50$, not $500$.

Here is my code:


tidsSteg = 100
kappaVektor = vector("numeric", tidsSteg)
sigmaVektor = vector("numeric", tidsSteg)
muVektor = vector("numeric", tidsSteg)
nrOfPoints = vector("numeric", tidsSteg)

for (i in 1:tidsSteg) {

  thomasSim = rThomas(kappa = 20, scale = 0.1, mu = 6, 
                      win = owin(c(0,1),c(0,1)),
                      nsim = 1, 
                      drop = TRUE, 
                      saveLambda = FALSE,
                      poisthresh = 1e-6, 
                      saveparents = FALSE)

  nrOfPoints[i] = thomasSim[["n"]]

  resultat = thomas.estK(thomasSim, startpar=c(kappa=20 ,scale=0.1), 
                         q = 1/4, p = 2, rmin = NULL, rmax = NULL)

  kappaVektor[i] = resultat$modelpar[[1]]
  sigmaVektor[i] = resultat$modelpar[[2]]
  muVektor[i] = nrOfPoints[i]/kappaVektor[i]
medelSigma = mean(sigmaVektor)
medelKappa = mean(kappaVektor)
medelMu = mean(muVektor)




Browse other questions tagged or ask your own question.