Как выровнять цветовые шкалы двух тепловых карт (функция image2D из пакета plot3D)?

Question

РЕДАКТИРОВАТЬ: Оказалось, что это «просто» неправильно поставленные скобки в densite_cauchy_norm, которые исказили результаты, теперь все работает нормально.
У меня есть три 2D гистограммы (тепловые карты), которые мне нужно сравнить.
По какой-то причине один из них (второй на прилагаемом изображении) отображает значения в цвета иначе, чем два других (например, плотность 0,020 на второй гистограмме имеет тот же цвет, что и плотность 0,015 на двух других).

Они получены при моделировании по методу Монте-Карло, где мне нужно сравнить гистограммы 2 выборок (X_1, X_2), смоделированных при отбраковке, с двумя огибающими различной плотности (графики 1 и 2).) к фактической «плотности» (f_tilde, ненормализовано), из которой они отбираются (график 3).Таким образом, пока все они не имеют одинакового цветового сопоставления, я не могу действительно показать, насколько они близки к теоретическому.

Для первых двух гистограмм я использовал одинаковые команды для обоих образцов X_1 и X_2.

Посмотрите на код ниже;Вы можете найти данные внизу:

library(plot3D)  
par(mfcol = c(1, 3)); 

n <- 100000;  
X_1 <- simu_f_1(n);  
X_2<- simu_f_2(n);  
x_s <- 50;  y_s <- 50;    

mon_histo(X_1, x_s, y_s, opt_3D = FALSE);  
mon_histo(X_2, x_s, y_s, opt_3D = FALSE);  

z <- outer(seq.int(0, 4, length= x_s), seq.int(0, 2, length = y_s), f_tilde);

image2D(z=z/sum(z), x = seq.int(0, 4, length= x_s),y = seq.int(0, 2, 
                    length = y_s), border="black", contour = FALSE);

Вторая гистограмма должна приближаться к наилучшей, но без той же цветовой шкалы, что и у других, она кажется неправильной.

Данные и необходимые функции для создания воспроизводимого примера:

densite_unif_norm <- function(x,y) { return(dnorm(x, mean = 2, sd = 1) * dunif(y, min = 0, max = 2)); }  
densite_cauchy_norm <- function(x,y) { return( dnorm(x, mean = 2, sd = 1)) * dcauchy(y, location = 1, scale = 0.5); }  
f_tilde <- function(x,y) {  
    return( (x>=0)*(x<=4)*(y>=0)*(y<=2) * exp(-0.5*(x-2)**2)*  
                (cos(x)**2+(2*sin(y)**2)*(cos(x)**4)) / (1+4*(y-1)**2)  );  
}
simu_f_1 <- function(n) { #simulates size n sample X_1 ~ "density" f_tilde by rejection sampling  
    M <- 6*sqrt(2*pi); # M = 15.04  
    f <- function(x) { 
        repeat{  
            Y <- c(rnorm(n = 1, mean = 2, sd = 1), runif(n = 1, min = 0, max = 2));  
            U <- runif(n = 1, min = 0, max = M*densite_unif_norm(Y[1], Y[2]));  
            if(U < f_tilde(Y[1], Y[2])) { break; } #la boucle s'arrête à la dernière valeur de y t.q u < f_tilde(y)  
        }  
        return(rbind(Y[1], Y[2]));  
    }
    return(Vectorize(f)(1:n));  
}

simu_f_2 <- function(n) { 
    M <- 3*(pi^1.5)/sqrt(2); # M = 11.81  
    f <- function(x) {  
        repeat{  
            Y <- c(rnorm(n = 1, mean = 2, sd = 1), rcauchy(n = 1, location = 1, scale = 0.5));  
            U <- runif(n = 1, min = 0, max = M*densite_cauchy_norm(Y[1], Y[2]));  
            if(U < f_tilde(Y[1], Y[2])) { break; }
        }  
        return(rbind(Y[1], Y[2]));  
    }  
    return(Vectorize(f)(1:n));  
}  
mon_histo <- function(X, x_sub, y_sub, opt_3D = FALSE) {  
    x_c <- cut(X[1,], x_sub);  
    y_c <- cut(X[2,], y_sub);  
    z <- table(x_c, y_c);  
    if(opt_3D) {hist3D(z=z, border="black"); }  
    else {   
        u <- seq(from = range(X[1,])[1], to = range(X[1,])[2], length.out = x_sub);  
        v = seq(from = range(X[2,])[1], to = range(X[2,])[2], length.out = y_sub);  
        image2D( z=z/length(X[1,]), border="black", x = u, y = v);  
    }  
}