kryl писал(а):Более строго можно переформулировать так:
Имеется область пространства от 0 до Xmax и от 0 до Ymax и общей площадью Smax=Xmax*Ymax (лесной фонд).
В ней имеется множество С1 С2 .. Сn кругов диаметром D общей площадью Sz случайно распределенных в пространстве(куртины усыхания). Внутри кругов значение усыхания Z=1, вне их – 0.
Имеется P1 P2 .. Pm прямоугольников(пробных площадей) размером a на b и общей площадью Spp случайно распределенных в пространстве.
Пересечение множества кругов С() и прямоугольников P() дает площадь усыхания внутри прямоугольников(пробных площадей) Szpp
Мы пытаемся оценить отношение Sz/Smax (площадь усохших насаждений в лесном фонде) зная соотношение Szpp/Spp (долю усохших деревьев на пробных площадях).
Sz/Smax=Szpp/Spp +/- Err (допустим при p=95%).
Задача найти функцию зависимости Err(ошибки) от D(диаметра куртин усыхания), Sz/Smax(доли усохших насаждений), a и b (размеров пробных площадей) и Spp/Smax(объема выборки).
В первом приближении можно заменить круги(куртины) квадратами.
тогда, если выборка достаточно большая, можно использовать бутстрап.
1) считаем выборку представительной, т.е. E[Sz/Smax] = E[Szpp/Spp], поэтому Szpp/Spp есть оценка
2) для оценки дисперсии делаем 10000 выборок с повтором из Szpp/Spp и вычисляем выборочное распределение средних, из него берем 5% границы
Код: Выделить всё
# 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))
Если мало - придется считать площади, а там еще больше мороки, сейчас на это времени нет (будет только в январе)
P.S. все это писалось "из головы", без проверки и отладки. Т.е. следует воспринимать как идею, а не работающую программу
