How to fix error in R code because of NaN values?

Question

I have an assignment to estimate parameter $θ$ from a sample with Pareto distribution with density $f(x; θ) = θ/x^(θ + 1), x ≥ 1$, where $θ>0$ is an unknown parameter. However, we do not know the realized sample $x$, we only know for each $x_i$ the given interval $(u_i; v_i)$ in which it is located.

Using the EM algorithm we have to estimate parameter $θ$. Also, EM has to be implemented in R and the code for that is down below.

When I run the code, I have an error because of NaN values. I've tried changing the starting value of the parameter, but NaN values still appear. How to fix this?

set.seed(1)

library(VGAM)
library(ggplot2)


#--------------------------------------------------------------------
# EM algorithm

# step E 
expected_xs <- function(teta, u, v) {
  teta/(teta-1) * 
    1/(1/(u^teta)-1/(v^teta))
}

# step M 
maximize_logL <- function(xs) {
  length(xs)/
    sum(log(xs))
}

EM_estimate <- function(teta_0, u, v, tol = 1e-8, maxiter = 1000) {
  xs <- expected_xs(teta_0, u, v)
  teta <- maximize_logL(xs)
  print(teta)
  iter <- 1
  
  while(!is.na(teta) && (abs(teta - teta_0) > tol) &&
        iter < maxiter) {
    iter <- iter + 1
    teta_0 <- teta
    
    xs <- expected_xs(teta_0, u, v)
    teta <- maximize_logL(xs)
    print(teta)
  }
  return(teta)
}
#--------------------------------------------------------------------




# Data

df <- read.table(header=T, text="
interval  freq
1 1-1.5    15
2 1.5-2    5
3 2-2.5    3 
4 2.5-3    3
5 3-1000   4") 
df


#u <- c(1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5)
#v <- c(1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10)
u <- seq(1, 999.5, by=0.5)
v <- seq(1.5, 1000, by=0.5)



teta1=EM_estimate(0.3, u, v)
teta1 

# we compare barplot with density (with its now estimated parameter)
barplot(df$freq, names.arg = df$interval)
curve(100*dpareto(x,teta1,1), add=TRUE, col="steelblue", lwd = 2)

One more thing, when I change teta/(teta-1)to teta/(teta+1) in here:

expected_xs <- function(teta, u, v) {
  teta/(teta-1) * 
    1/(1/(u^teta)-1/(v^teta))
}

everything works normally.