# # Función de detección # (probabilidad de detectar un objeto a distancia "d" del transecto # g=function(x) { g=1-0.0052*x^2 g=ifelse(g<0,0,g) return(g) } # # Función para decidir si un objeto a distancia "d" ha sido detectado o no. # Devuelve 1 si el objeto se detecta, 0 si no se detecta. # Puede aplicarse a un vector "d" de distancias # detect=function(d) { det01=function(x) {rbinom(1,1,g(x))} return(sapply(d,det01)) } # # Dibujamos la función de detección entre 0 y k km: k=25 base=seq(0,k,by=0.05) gb=g(base) plot(base,gb,type="l",ylim=c(0,1),col="red") # # El área a investigar es un cuadrado de área l^2 km2 con centro en (0,0): # l=100 # # D es la densidad de objetos por km2 # D=0.4 # # Distancia máxima de detección # delta=40 # # El número de objetos se simula de acuerdo a una Poisson de tasa Dx(l^2) # lambda=D*(l^2) N=rpois(1,lambda) # # Se generan las posiciones de los N objetos: # # El centro es (0,0), y el cuadrado va de -l/2 a l/2 x=runif(N,-l/2,l/2) y=runif(N,-l/2,l/2) objetos=data.frame(cbind(x,y)) plot(objetos$x,objetos$y) # # Trazamos un transecto aleatorio # Para ello elegimos un ángulo entre 0 y 2*pi y trazamos un segmento con esa # pendiente asociada; la longitud del segmento es l/4 # phi=runif(1,0,2*pi) p1=(l/4)*c(cos(phi),sin(phi)) p2=(l/4)*c(cos(phi+pi),sin(phi+pi)) extremos=rbind(p1,p2) # puntos extremos del transecto lines(extremos,col="red") # se dibuja el transecto # # A continuación se seleccionan los puntos que se encuentran a una distancia # inferior a delta del transecto # ord1=p1[2]+(cos(phi)/sin(phi))*p1[1] ord2=p2[2]+(cos(phi)/sin(phi))*p2[1] ppt=-(cos(phi)/sin(phi)) # pendiente recta perpendicular al transecto abline(ord1,ppt,col="green") # estas rectas definen los límites de abline(ord2,ppt,col="green") # observación perpendiculares al transecto # dist=abs(x*sin(phi)-y*cos(phi)) # Distancia de cada punto a la recta # Los puntos con s<=0 son los que están entre los extremos del transecto s=sign((y-p1[2]-ppt*(x-p1[1]))*(y-p2[2]-ppt*(x-p2[1]))) objetos=cbind(objetos,dist,s) # # Pintamos de rojo los puntos detectados: cerca=subset(objetos,((dist