library(spatstat) ########################## ### Minimalni kontrast ### ########################## # Modely shlukovych bodovych procesu muzeme fitovat na data metodou minimalniho kontrastu # pomoci prikazu "kppm". Defaultne se vyuziva K-funkce. # Nafitujme model Thomasova procesu na dataset "redwood" (shlukovy proces): data(redwood) fit <- kppm(redwood, ~1, "Thomas") fit plot(fit) # Prvni argument funkce kppm je zkoumany bodovy vzorek; # druhy argument je vzorec popisujici logaritmus funkce intenzity v modelu (~1 odpovida stacionarite); # treti argument udava typ shlukovani (ruzne typy popsane slovne). # Nafitujme model Maternova shlukoveho procesu: fitM <- kppm(redwood, ~1, "MatClust") fitM plot(fitM) # Vysledek je objekt tridy "kppm". # Z nafitovaneho modelu muzeme rovnou simulovat (diky funkci simulate.kppm): plot(simulate(fit, nsim = 4)) # Muzeme take vyrabet obalky odpovidajici odhadnutemu modelu pomoci envelope.kppm: plot(envelope(fit, Lest, nsim = 39)) # Metodu minimalniho kontrastu muzeme vyuzit i pro jine charakteristiky, napr. pro parovou # korelacni funkci (vic ale neni ve Spatstatu implementovano): fitp <- kppm(redwood, ~1, "Thomas", statistic = "pcf") fitp attributes(fitp) # Podrobnejsi informace napr. o nastaveni odhadovaci procedury (integracni obor, exponenty p, q): fitp$Fit # Tyto parametry muzeme nastavovat sami, napr. na zaklade nejakych doporuceni z literatury: fitp2 <- kppm(redwood, ~1, "Thomas", statistic = "pcf", q=1/2) fitp2 fitp2$Fit fitp3 <- kppm(redwood, ~1, "Thomas", statistic = "pcf", q=1/2, rmax=0.15) fitp3 fitp3$Fit # Metodu min. kontrastu muzeme vyuzit i pro jine modely, pro nez zname teoreticky predpis # pro K-funkci ci parovou korelacni funkci, jako funkce nejakych parametru, # treba alespon ve forme integralniho vzorce. # Prikladem je log-gaussovsky Coxuv proces (s exponencialni kovariancni funkci): ?lgcp.estK # min. kontrast na K-funkci ?lgcp.estpcf # min. kontrast na parovou korelacni funkci # V tomto pripade musime inicializovat optimalizacni proceduru (zadat pocatecni odhad parametru): fitLGCP <- lgcp.estK(redwood, c(sigma2 = 0.1, alpha = 1)) fitLGCP fitLGCP2 <- lgcp.estpcf(redwood, c(sigma2 = 0.1, alpha = 1), rmin=0.01, rmax=0.25) fitLGCP2 ########################### ### Slozena verohodnost ### ########################### # Pro shlukove procesy muzeme pouzit take metodu slozene verohodnosti: fitCL <- kppm(redwood, ~1, "Thomas", method="clik") fitCL ######################### ### Pseudoverohodnost ### ######################### # Uvazujeme konecny bodovy proces s hustotou, na omezenem okne, ktery ma "rozumny tvar" # Papangelouovy podminene intenzity. # K odhadu parametru pouzijeme prikaz "ppm" ve tvaru "ppm(X, ~ terms, V)". # Prvnim argumentem je bodovy vzorek; # druhy argument popisuje cleny prvniho radu v potencialu; # treti argument je objekt tridy "interact" a popisuje typ interakci druheho radu. # Nafitujme Straussuv proces na "cells" dataset (prepdokladame, ze interakcni vzdalenost r zname): data(cells) fitCELLS <- ppm(cells, ~1, Strauss(r = 0.1)) fitCELLS # Muzeme uvazovat nehomogenitu ve clenech prvniho radu, napr. log-linearni zavislost na souradnicich: fitCELLS2 <- ppm(cells, ~x + y, Strauss(r = 0.1)) fitCELLS2 # Jsou k dispozici i jine modely, napr. Geyeruv saturacni proces (jde pouzit i k modelovani shluku): ppm(redwood, ~1, Geyer(r = 0.07, sat = 2)) # Z nafitovaneho modelu muzeme rovnou simulovat pomoci prikazu "rmh" (vyuziva MCMC retezec): fitCELLS <- ppm(cells, ~1, Strauss(r = 0.1)) Xsim <- rmh(fitCELLS) par(mfrow=c(1,2), mar=c(1,1,1,1)) plot(cells, main = "Cells") plot(Xsim, main = "Simulation") dev.off() # zavre graficky vystup, soucasne vrati defaultni nastaveni pro pristi obrazek. # Kresleni obalek: plot(envelope(fitCELLS, nsim = 39)) # Rusive parametry, jako napr. dosah interakci "r" u Straussova procesu, nemohou byt primo odhadnuty # pomoci prikazu "ppm". Ani teorie k odhadu techto parametru neni hotova (uspokojiva). # V nekterych specialnich pripadech je znam maximalne verohodny odhad rusiveho parametru, napr. # pro "hard-core proces" (Straussuv proces s interakcnim parametrem gama = 0) je rusivym parametrem # dosah interakci "r" a jeho MLE odhadem je nejmensi pozorovana vzdalenost mezi dvema body. # Ma tedy smysl napriklad pro "cells" dataset postupovat takto: rhat <- min(nndist(cells)) rhat <- rhat * 0.99999 ppm(cells, ~1, Strauss(r = rhat)) # Muzeme ale vyuzit tzv. profile pseudolikelihood (hledani rusiveho parametru na nejake siti), kdy # pro kazdou uvazovanou hodnotu rusiveho parametru maximalizujeme pseudoverohodnost pres ostatni, # regularni parametry. data(simdat) # simulovany dataset df <- data.frame(r = seq(0.05, 2, by = 0.025)) # vybirame uvazovane hodnoty "r" df pfit <- profilepl(df, Strauss, simdat, ~1) # hlavni vypocet # Shrnuti: pfit plot(pfit) # Nejlepsi dostupny odhad parametru: pfit$fit # Ukol: na dataset "swedishpines" nafitujte vhodny model Straussova procesu a ukazte # pomoci obalkoveho a odchylkoveho testu (dclf.test, mad.test), ze data odpovidaji nafitovanemu modelu. data(swedishpines) plot(swedishpines) ?swedishpines