Форумы GIS-Lab.info

Код: Выделить всё

x<-seq(0,100,by=0.1)
y<-2*sqrt(50^2-(x-50)^2)
y<-y/sum(y)
p<-expand.grid(x1=1:length(x),x2=1:length(x))
p.g<-matrix(y[p[,1]]*y[p[,2]],ncol=length(x),nrow=length(x))
p.range<-range(p.g)
image(x, x, p.g, col = terrain.colors(100), axes = FALSE)
contour(x, x, p.g, levels = seq(p.range[1], p.range[2], length.out = 10), add = TRUE, col = "peru")

Код: Выделить всё

# Szpp_Spp.dat - площадь пересечения одной куртины с одной  пробной площадкой, число строк равно числу проб
# нули тоже пишем в данные, заголовка нет
p<-drop(read.table("Szpp_Spp.dat",header=FALSE))
n.data<-length(p)
n.boot<-10000
p.stat<-rep(NA,n.boot)
for(i in 1:n.boot) {
  ind<-sample(1:n.data,replace=TRUE)
  p.stat[i]<-mean(p[ind])
}
print(quantile(p.stat,c(0.025,0.975)))

Код: Выделить всё

p.D <-20 # "диаметр" куртины (куртина - это квадрат)
p.lx<-5 # размер пробной площадки по X
p.ly<-3 # размер пробной площадки по Y
p.Xmax<-500
p.Ymax<-800

# считаем без учета всяких тонкостей. Считаем p.D << min(p.Xmax,p.Ymax), попасть в 2 курины невозможно

P.x<-(p.D + p.lx)/(p.Xmax-p.lx) # вероятность "зацепить" пробной площадкой прекцию куртины по X
P.y<-(p.D + p.ly)/(p.Ymax-p.ly) # вероятность "зацепить" пробной площадкой прекцию куртины по Y
P.xy<-P.x * P.y # вероятность "зацепить" куртину пробной площадкой

p.fact<-as.integer(p > 0) # какими площадками зацепили
p.y<-sum(p.fact)          # сколько площадок зацепились за куртины

# логарифмическая функция правдоподобия (вероятность зацепить одну из k куртин есть k*P.xy в силу редкости куртин):
# LL(p.y)=p.y*ln(p.k*P.xy) + (n.data-p.y)*ln(1-p.k*P.xy)
# находим максимум по k
k.ML<-p.y/(P.xy*n.data)
S.ML<-k.ML*p.D^2 / (p.Xmax*p.Ymax)

# дисперсия числа куртин K: Var(K)=n.data*(1-K*P.xy)*(K*P.xy)
k.ML.var<-n.data*(1-k.ML*P.xy)*(k.ML*P.xy)

# дисперсия S.ML: Var(S.ML)=Var(k)*(p.D^2 / (p.Xmax*p.Ymax))^2
S.ML.var<-k.ML.var*(p.D^2 / (p.Xmax*p.Ymax))^2

# считая биномиальное распределение более-менее нормальным, получаем интервал:
S.ML.sd<-sqrt(S.ML.var)
print(c(S.ML+qnorm(0.025,0,1)*S.ML.sd,S.ML+qnorm(0.975,0,1)*S.ML.sd))

Код: Выделить всё

k<-100
rez1<-matrix(0,nrow=k,ncol=1)
rez2<-matrix(0,nrow=k,ncol=1)

for (k in 1:100){
Xmax<-800
us<-matrix(0,nrow=Xmax,ncol=Xmax)
pp<-matrix(0,nrow=Xmax,ncol=Xmax)
n<-200
m<-200
D<-10
plx<-5
ply<-3
Sz<-D^2*n
Spp<-plx*ply*m
Smax<-Xmax^2
Dz<-(Sz/Smax)
Dzpp<-(Szpp/Spp)

for (i in 1:m) {
x_rand<- ceiling(runif(1)*Xmax)
y_rand<-ceiling(runif(1)*Xmax)
   for(i3 in 1:ply) {
    for(i2 in 1:plx) {
	xi<-x_rand+i2-1
	yi<-y_rand+i3-1
	xic=ifelse(xi<Xmax, xi, xi-Xmax)
        yic=ifelse(yi<Xmax, yi, yi-Xmax)
    pp[(xic), (yic)]<-1
    }
   }
}
for (i in 1:n) {
x_rand<- ceiling(runif(1)*Xmax)
y_rand<-ceiling(runif(1)*Xmax)
   for(i3 in 1:D) {
    for(i2 in 1:D) {
	xi<-x_rand+i2-1
	yi<-y_rand+i3-1
	xic=ifelse(xi<Xmax, xi, xi-Xmax)
        yic=ifelse(yi<Xmax, yi, yi-Xmax)
    us[(xic), (yic)]<-1
    }
  }
}
inter1<- us*pp
Szpp<-sum(inter1)
Err<-(Szpp/Spp)-(Sz/Smax)


rez1[k,1]<-Szpp
rez2[k,1]<-Err
}

Код: Выделить всё

p.fact<-as.integer(p > 0) # какими площадками зацепили