需要积分知识

等于pi*r0^2*D0=29.45?


这样就算出来了，还要模拟干什么？

我的研究有些复杂。比如样地的昆虫分布密度在样地是均匀的。坐标轴的样点为陷阱（含有引诱物），用以捕捉昆虫。但是相邻的陷阱会对周围昆虫产生影响，比如介于两个引诱半径之间，可能去一个也可能去另外一个。

在实际中希望得到昆虫的种群密度，而陷阱引诱半径是不知道的，于是要测算两个参数。引诱的半径和在此引诱圆圈中昆虫的点密度。

具体就不跟你讲，只需要看我的问题就可以了，谢谢！

帮顶

我的问题，其实很简单。一个正方形，假设点的数量是一定的，随机分布在这个正方形中。在正方形中有一个小圆，希望小圆内的点的密度正好等于（或者接近）= 正方形内点的总数量除以正方形的面积（即正方形中点的密度）。

大家可以使用R中的随机函数试试看，结果总是不等的。问题解决了，其实是相等的。呵呵。


# Form the point array of traps in the right x-axis
#x.right <- c()
#temp    <- 0.5
#for (i in 1:12){
#    x.right <- c(x.right, temp)
#    temp    <- temp + i
#}

x.right <- c()
temp    <- 0.5
for (i in 1:7){
    x.right <- c(x.right, temp)
    temp    <- temp + 2^(i-1)
}     

     
#x.right <- x.right[-1]


x.left  <- -x.right
y.left  <- rep(0, length(x.left))
y.right <- rep(0, length(x.right))

y.up    <- x.right
y.down  <- x.left
x.up    <- rep(0, length(y.up))
x.down  <- rep(0, length(y.down))

# x and y are used to store the planar coordinates of all traps
x       <- c(x.left, x.right, x.up, x.down, 0)
y       <- c(y.left, y.right, y.up, y.down, 0)




# Presume that the density and radius are known
D0 <- 1.5
r0 <- 2.5




# Give a rectangle of w * w
w  <- max(x.right) + r0


# 找到错误了，模拟的面积错误！应该为 A1 <- ( 2 * w ) * ( 2 * r0 ); A2 <- ( 2 * r0 ) * ( 2 * w )
A1 <- w*(2 * r0)
A2 <- (2 * r0)*w
##########################################################


A3 <- (2 * r0)^2

A <- (A1 + A2 - A3)
point.number <- A * D0
#point.number

q <- 7

n.sim       <- 60
distance    <- c()
number.circled <- c()

for(h in 1:n.sim){
     cat(" ", h)
     if (h %% 10 == 0)
     cat("\n")





     # Insect points are randomly distributed on the areas around the axes
     x1.point <- runif(point.number/2, -w, w)
     y1.point <- runif(point.number/2, -r0, r0)
     x2.point <- runif(point.number/4, -r0, r0)
     y2.point <- runif(point.number/4, -w, -r0)
     x3.point <- runif(point.number/4, -r0, r0)
     y3.point <- runif(point.number/4, r0, w)


     x.point <- c(x1.point, x2.point, x3.point)
     y.point <- c(y1.point, y2.point, y3.point)

     #x.point <- runif(w^2*D0, -w, w)
     #y.point <- runif(w^2*D0, -w, w)



     sign.array <- c()
     x.array    <- c()
     y.array    <- c()

     # The first circle is for all the small trap circles with radius r0  
     for ( j in q ){      
         # The second circle is for all the insect points
         for ( i in 1:length(x.point) ){
             if (x.point>= x[j] - r0 & x.point <= x[j] + r0 & y.point >= y[j] - sin(acos(abs(x.point
                        -x[j])/r0))*r0 & y.point <= y[j] + sin(acos(abs(x.point-x[j])/r0))*r0){
                 sign.array <- c(sign.array, j)
                 x.array    <- c(x.array, x.point)               
                 y.array    <- c(y.array, y.point)
             }               

         }
     }
     
number.circled <- c(number.circled, length(sign.array))

}

mean(number.circled)
sd(number.circled)




## An example of trap circle with insect points included
x.range <- c()
y.range <- c()
step    <- 2*pi / 1000
theta   <- seq(0, 2*pi - step, step)
for (i in 1:length(theta)){
    x.range <- x[q] + cos(theta)*r0
    y.range <- y[q] + sin(theta)*r0
}



plot(x.array, y.array, xlim=c(x[q]-r0, x[q]+r0), ylim=c(y[q]-r0, y[q]+r0), xlab=c("x"), ylab=c("y"), cex = 2)
points(x[q], y[q], pch=16, col=2)
points(x.point, y.point, pch=4, cex = 2 )
lines(x.range, y.range)
length(x.array)

x11()
plot(x.point, y.point, xlim =c(x[q]-r0, x[q]+r0), ylim = c(y[q]-r0, y[q]+r0), xlab=c("x"), ylab=c("y"), cex=2)