Optimization with a function that involves integration

Question

I'm working on a problem where a parameter is estimated through minimizing the sum of squares. The equations needed are:

I used optim in the package stats:

# provide the values for a test dataset (the y estimated should be 1.41)
pvector <- c(0.0036,0.0156,0.0204,0.0325,0.1096,0.1446,0.1843,0.4518)
zobs <- c(0.0971,0.0914,0.1629,0.1623,0.3840,0.5155,0.3648,0.6639)
# make input of the C value 
c <- function(y){
  gamma(y)/((gamma(y*(1-pvector)))*(gamma(y*pvector)))
}
# make input of the gamma function
F1 <- function(y){
  f1 <- function(x){
    c*(1-x)^(y*(1-pvector)-1)*x^(y*pvector-1)
  }
  return (f1)
}
# integration over x
int <- function(y){
  integrate (F1(y),lower =0.001, upper =1) 
}
# write the function for minimization 
f2 <- function(y) {
  sum ((int-zobs)^2)
}
# minimization 
optim(0.01,f2, method = "Brent", lower =0, upper = 1000, hessian=TRUE)

Which didn't work. I received the following error message:

Error in int - zobs : non-numeric argument to binary operator

I think there must be something fundamentally wrong with the way how the function was written.