how to write the joint pmf of bivariate distribution which defined over more than one interval in r

Question

definition of the bivariate distribution If I have the following joint probability density function: The joint PMF i.e, the joint probability mass function (pmf) is consisting of the following pmf and cumulative distribution function The marginal pmf

The R code of the marginal distribution is given as follows (i.e the pmf of DIW distribution)

ddiw<- function(x, t, eta){  # t means theta parameter
  stopifnot( eta>0, x>=0)
  pmf<-t^(x+1)^(-eta) - t^(x)^(-eta)
  return(pmf)
}

Its cumulative distribution function in R is as follows

pdiw<-function(x, t, eta){
  stopifnot(  eta > 0)
  cdf<- t^(x+1)^(-eta)
  return(cdf) 
}

I want to write the joint pmf in R as in equation (4). I tried to write the joint pmf in equation (4) in R, but I did not succeed. Also, I want to plot the joint pmf in R as in the following figure 3D plot of bivariate distribution

could you help me to write the joint pmf in R and plot it as in the given figure.

Thanks in advance.

Edit

I write the joint pmf to be more clear as follows The joint pmf

The joint pmf when substituting the pmf of DIW and substituting cdf of DIW

Answer 1

Not sure if this is helpful, but you could write the function along the lines:

joint_pmf <- function(x1, x2, eta) {
    stopifnot(
        all(x1>0), all(x2>0), 
        all(is.finite(x1)), all(is.finite(x2)),
        length(x1)==length(x2)
    )
    # the result vector
    n <- length(x1)
    pmf <- rep(NA, n) 

    # check for the first condition
    idx <- which(x1<x2)
    pmf[idx] <- NA # here you need to fill in

    # check for the second condition
    idx <- which(x1>x2)
    pmf[idx] <- NA # here you need to fill in

    # and so on..
}

This function accepts vectors for the x1, x2 arguments. Since the pmf is defined piecewise, the result is computed for subsets only.