Fitting of bivariate distribution using R

Question

#Bivariate data
#########LIPD##############
library(fitdistrplus)
library(AdequacyModel)
library(extRemes)
library(hypergeo)
library(gsl)
library(fitdistrplus)
#
x<-c(1:53)
#
y<-c(0.054, 0.064, 0.704, 0.816, 0.235, 0.976, 0.865, 0.364, 0.479, 0.568, 0.352, 0.978, 0.787, 0.976, 0.087, 0.548, 0.796, 0.458,
     0.087, 0.437, 0.421, 1.978, 1.756, 2.089, 2.643, 2.869, 3.867, 3.890, 3.543, 3.079, 3.646, 3.348,
     4.093, 4.092, 4.190, 4.237, 5.028, 5.083, 6.174, 6.743, 7.274, 7.058, 8.273, 9.324, 10.827, 11.282,
     13.324, 14.278, 15.287, 16.978, 17.209, 19.092, 20.083)
#

dlip <- function(x, y ,l,a,r){
  
  (l * (l+(r*x))^(x-1) * y^(x+1) * exp(-((r*x*y)+(y/a))))/(a^2 * factorial(x))
}
#
plip <- function(q,m,l,a,r){
  cumsum(dlip(q,m,l,a,r))
}
#
plip(x, y, 2,0.2,0.3)
#
est1 <- fitdist(x, y, "lip", start = list(l = 0.1,a = 0.1, r = 0.1),
                optim.method = "Nelder-Mead", lower = c(0, 0, 0), upper = c(Inf, Inf, 1))
est1
summary(est1)
gofstat(est1)
#

I have tried the above R code for finding mle of my bivariate mixture distribution.

In the R code l, a, and r are the parameters of the my bivariate distribution. But my R code cannot be run.

Can anyone help me to solve the problem?