准备画次序统计量的二维联合密度函数 在最后绘图的时候 程序报错： 被替换的项目不是替换值长度的倍数
了解到是因为函数中有f[x-y>=0]，f[x-y<0]这么一个需要比较大小的步骤才出现这个报错
但是不知道具体该怎么改进
希望各位能给出指导

感激不尽！
程序如下：
n=1000
y=rnorm(n,0,1)
y1=sort(y)
k<-20
j<-300  #选择k=20，j=300
fun.cixu2=function(x,y)
{A=c(n,k-1,j-k-1,n-j)
B=c(0,0,0,0)
for (l in 1:4) {
  for (i in 1:A[l]) {
    B[l]=B[l]+log(i)
  }}
b=exp(B[1]-B[2]-B[3]-B[4]) #通过转化计算常数项的值            
f=rep(NA,length(x))
f[x-y<0]=b*pnorm(x,0,1)^(k-1)*(1-pnorm(x,0,1))^(n-j)*dnorm(x,0,1)*dnorm(y,0,1)*(pnorm(y,0,1)-pnorm(x,0,1))^(j-1-k) #判断x，y的大小并代入计算
f[x-y>=0]=0 #若x>=y，则函数值为0
return(f)
}
fun.cixu2(y1[20],y1[300]) #取之前得到的随机数进行代入计算
x <- seq(-3,1,length.out = 100)
y <- seq(-3,1,length.out = 100)
z= outer(x,y,fun.cixu2)
persp(x=seq(-3,1,length.out = 100),y=seq(-3,1,length.out = 100),z,theta=-25,phi=25,col='lightblue') #画出二维联合密度图像