一个奇怪的模拟问题，请赐教

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

现在我遇到这么一个头疼的问题：

在平面坐标系中，假设沿着x正方向依次有0.5, 1.5, 3.5, 7.5, 15.5, 31.5, 63.5. 现在假设以最后一个点为圆点，以r0 = 2.5为半径做一个圆。
假设沿着x轴正方向有一个长70，宽2倍r0的矩形。假设这个矩形内有若干点，点的密度为D0 = 1.5，那么点的数量其实就是D0*70*2*r0（也即，密度乘以面积A=70*2*r0）。
那么请问落在上述圆内的点的数量是多少呢？是不是等于pi*r0^2*D0=29.45?
但是我怎么模拟都是14.75左右。我非常奇怪。我使用了R和Matlab软件结果都一样。下边是Matlab程序。请高手指点迷津，谢谢！

w  = 70;
r0 = 2.5;
D0 = 1.5;

% 找到错误了，下述代码应该改为： A = (2*w)*(2*r0);
A = w*(2*r0);

point_number = A*D0;
temp = 0.5;
for i=1:7
    x_right(i) = temp;
    temp = temp + 2^(i-1);
end
x = x_right;
y = zeros(length(x), 1)';
n_sim = 200;

q = 1;
for i=1:n_sim
   
    i
    x_point = - w + (2 * w) * rand(point_number, 1);
    y_point = - r0 + (2 * r0) * rand(point_number, 1);
    %plot(x_point, y_point,'ko')
    counter = 0;
    for j=1:length(x_point)
        if x_point(j) >= x(q) - r0 & x_point(j) <= x(q) + r0 & y_point(j) >= y(q) - sin(acos(abs(x_point(j) - x(q))/r0))*r0 & y_point(j) <= y(q) + sin(acos(abs(x_point(j) - x(q))/r0))*r0
            counter = counter + 1;
        end
    end
    Number(i) = counter;
end
figure(1)
plot(x_point, y_point,'ko')
axis square
axis([x(q)-r0 x(q)+r0 y(q)-r0 y(q)+r0])
mean(Number)

k = 0;
s = 10000;
for theta = 0:(2*pi/s):(2*pi- 2*pi/s)
    k = k + 1;
    x_temp(k) = x(q) + cos(theta)*r0;
    y_temp(k) = y(q) + sin(theta)*r0;
end
%%%% An example %%%%
figure(2)
plot(x_point, y_point,'ko')
hold on
plot(x(q), y(q), 'r+')
plot(x_temp, y_temp, '-')
axis([x(q)-r0 x(q)+r0 y(q)-r0 y(q)+r0])
axis square
hold off

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：MATLAB程序 matlab软件 counter Number length 坐标系 程序 平面 软件 左右

需要积分知识

等于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)

两点供参考:
   1. x_point = - w + (2 * w) * rand(point_number, 1);
      范围取大了, -w ~ w

   2.acos(abs(x_point(j) - x(q))/r0)会产生complex
     For real elements of X in the domain[-1,1], acos(X) is real and in the range[0,pi]
     For real elements of X outside the domain , acos(X) is complex.
     换句话说x_point(j)要先取范围内的值

simulation=2000,可得mean(Number)=29.4955

呵呵，我早晨发现了原因，原来是我设定的矩形面积应该为A=(2w)*(2*r0)，点数为A*D0，而我则错误地只投射了w*(2*r0)个点即A*D0/2，所以得到的点数只有原来的一半。谢谢！