Last updated on 2020-02-19 10:49:01 CET.
| Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
|---|---|---|---|---|---|---|
| r-devel-linux-x86_64-debian-clang | 0.1.0 | 2.68 | 19.96 | 22.64 | ERROR | |
| r-devel-linux-x86_64-debian-gcc | 0.1.0 | 1.89 | 16.23 | 18.12 | ERROR | |
| r-devel-linux-x86_64-fedora-clang | 0.1.0 | 28.47 | ERROR | |||
| r-devel-linux-x86_64-fedora-gcc | 0.1.0 | 27.71 | ERROR | |||
| r-devel-windows-ix86+x86_64 | 0.1.0 | 6.00 | 33.00 | 39.00 | OK | |
| r-devel-windows-ix86+x86_64-gcc8 | 0.1.0 | 10.00 | 52.00 | 62.00 | OK | |
| r-patched-linux-x86_64 | 0.1.0 | 2.03 | 22.50 | 24.53 | OK | |
| r-patched-solaris-x86 | 0.1.0 | 43.60 | OK | |||
| r-release-linux-x86_64 | 0.1.0 | 1.85 | 22.06 | 23.91 | OK | |
| r-release-windows-ix86+x86_64 | 0.1.0 | 8.00 | 35.00 | 43.00 | OK | |
| r-release-osx-x86_64 | 0.1.0 | OK | ||||
| r-oldrel-windows-ix86+x86_64 | 0.1.0 | 3.00 | 34.00 | 37.00 | OK | |
| r-oldrel-osx-x86_64 | 0.1.0 | OK |
Version: 0.1.0
Check: examples
Result: ERROR
Running examples in 'PACBO-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: PACBO
> ### Title: PACBO
> ### Aliases: PACBO
>
> ### ** Examples
>
> ## generating 4 clusters of 100 points in \strong{R}^{5}.
> set.seed(100)
> Nb <- 4
> d <- 5
> T <- 100
> proportion = rep(1/Nb, Nb)
> Mean_vectors <- matrix(runif(d*Nb,min=-10, max=10),nrow=Nb,ncol=d, byrow=TRUE)
> mydata <- matrix(replicate(T, rmnorm(1, mean= Mean_vectors[sample(1:Nb, 1, prob = proportion),],
+ varcov = diag(1,d))), nrow = T, byrow=T)
> R <- max(sqrt(rowSums(mydata^2)))
> ##run the algorithm.
> result <- PACBO(mydata, R, plot_ind = TRUE)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
PACBO
--- call from context ---
PACBO(mydata, R, plot_ind = TRUE)
--- call from argument ---
if (class(mydata) != "matrix") {
mydata = matrix(mydata, nrow = 1)
}
--- R stacktrace ---
where 1: PACBO(mydata, R, plot_ind = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE TRUE
--- function from context ---
function (mydata, R, coeff = 2, K_max = 50, scaling = FALSE,
var_ind = FALSE, N_iterations = 500, plot_ind = FALSE, axis_ind = c(1,
2))
{
if (class(mydata) != "matrix") {
mydata = matrix(mydata, nrow = 1)
}
if (R < max(sqrt(rowSums(mydata^2)))) {
print(c("R should be bigger than the maximum Euclidean distance of observations"))
}
else {
if (scaling) {
mydata = scale(mydata)
}
d = length(mydata[1, ])
T = length(mydata[, 1])
multiplier_R = 1.5
Nclusters = rep(1, N_iterations)
Niter_centers = list()
parameter_means_proposal = list()
sum_loss = rep(0, N_iterations)
if (var_ind) {
nb_of_clusters = rep(1, T)
pred_centers = list()
predicted_loss = rep(0, T)
proposal_loss = rep(0, T)
lambda_1 = rep(1, T)
lambda_1[2:T] = sapply(2:T, function(t) 2.5 * coeff *
(d + 2) * (t - 1)^(-0.5)/R)
lambda_2 = sapply(1:T, function(t) (d + 2) * (t)^(-0.5)/(multiplier_R *
R)^4)
c_1 = runiform_ball(5000, d, multiplier_R * R)
index_1 = which(apply((c_1 - t(replicate(5000, mydata[1,
])))^2, 1, sum) == instantaneous_loss(c_1, mydata[1,
]))
pred_centers[[1]] = c_1[index_1, ]
predicted_loss[1] = instantaneous_loss(pred_centers[[1]],
mydata[1, ])
pb <- txtProgressBar(0, 1, char = "=")
for (t in 2:T) {
tau_proposal = (K_max * t * d)^(-0.5)
Nclusters[1] = nb_of_clusters[t - 1]
parameter_means_proposal[[1]] = kmeans(matrix(mydata[1:(t -
1), ], ncol = d), Nclusters[1], nstart = 2,
iter.max = 10)$centers
Niter_centers[[1]] = t(apply(parameter_means_proposal[[1]],
1, function(x) rmt(1, mean = x, S = diag(tau_proposal,
d), df = 3)))
while (sum(sqrt(rowSums((Niter_centers[[1]]^2))) >
multiplier_R * R) > 0) {
Niter_centers[[1]] = t(apply(parameter_means_proposal[[1]],
1, function(x) rmt(1, mean = x, S = diag(tau_proposal,
d), df = 3)))
}
proposal_loss[1:(t - 1)] = apply(matrix(mydata[1:(t -
1), ], nrow = t - 1), 1, function(x) instantaneous_loss(Niter_centers[[1]],
x))
sum_loss[1] = sum(proposal_loss[1:(t - 1)]) +
0.5 * sum(lambda_2 * (proposal_loss - predicted_loss)^2)
for (n in 2:N_iterations) {
proposal_loss_temp = rep(0, T)
transition_prob = c(transition_probability(Nclusters[n -
1], Nclusters[n - 1] - 1, K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1], K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1] + 1, K_max))
transition_prob = transition_prob/sum(transition_prob)
new_k = sample(c(Nclusters[n - 1] - 1, Nclusters[n -
1], Nclusters[n - 1] + 1), 1, prob = transition_prob)
if (new_k == Nclusters[n - 1]) {
m_t = parameter_means_proposal[[n - 1]]
}
else {
if (new_k >= t - 1) {
new_k = t - 1
m_t = matrix(mydata[1:(t - 1), ], ncol = d,
byrow = F)
}
else {
m_t = kmeans(matrix(mydata[1:(t - 1), ],
ncol = d), new_k, nstart = 2, iter.max = 10)$centers
}
}
c_k_prime = t(apply(m_t, 1, function(x) rmt(1,
mean = x, S = diag(tau_proposal, d), df = 3)))
if (sum(sqrt(rowSums(c_k_prime^2)) < multiplier_R *
R) == new_k) {
log_numerator_prop = log(apply(matrix(1:Nclusters[n -
1], ncol = 1), 1, function(x) dmt(Niter_centers[[n -
1]][x, ], mean = parameter_means_proposal[[n -
1]][x, ], S = diag(tau_proposal, d), df = 3)))
log_denominator_prop = log(apply(matrix(1:new_k,
ncol = 1), 1, function(x) dmt(c_k_prime[x,
], mean = m_t[x, ], S = diag(tau_proposal,
d), df = 3)))
log_numerator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), new_k)
log_denominator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), Nclusters[n - 1])
ln_division = sum(c(log_numerator_prior,
log_numerator_prop) - c(log_denominator_prior,
log_denominator_prop))
proposal_loss_temp[1:(t - 1)] = apply(matrix(mydata[1:(t -
1), ], nrow = t - 1), 1, function(x) instantaneous_loss(c_k_prime,
x))
s_loss_prime = sum(proposal_loss_temp) +
0.5 * sum(lambda_2 * (proposal_loss_temp -
predicted_loss)^2)
ln_accept_ratio = (-lambda_1[t - 1] * (s_loss_prime -
sum_loss[n - 1])) + ln_division + log(transition_probability(new_k,
Nclusters[n - 1], K_max)) - log(transition_probability(Nclusters[n -
1], new_k, K_max))
}
else {
ln_accept_ratio = log(0)
}
bool = (log(runif(1)) < ln_accept_ratio)
if (bool) {
Niter_centers[[n]] = c_k_prime
parameter_means_proposal[[n]] = m_t
Nclusters[n] = new_k
sum_loss[n] = s_loss_prime
}
else {
Niter_centers[[n]] = Niter_centers[[n - 1]]
parameter_means_proposal[[n]] = parameter_means_proposal[[n -
1]]
Nclusters[n] = Nclusters[n - 1]
sum_loss[n] = sum_loss[n - 1]
}
}
nb_of_clusters[t] = Nclusters[N_iterations]
pred_centers[[t]] = Niter_centers[[N_iterations]]
predicted_loss[t] = instantaneous_loss(pred_centers[[t]],
mydata[t, ])
setTxtProgressBar(pb, t/T)
}
labels = labels_function(pred_centers[[T]], mydata)
if (plot_ind) {
plot(mydata[, axis_ind[1]], mydata[, axis_ind[2]],
xlim = c(min(mydata[, axis_ind[1]]) - 5, max(mydata[,
axis_ind[1]]) + 5), ylim = c(min(mydata[,
axis_ind[2]]) - 5, max(mydata[, axis_ind[2]]) +
5), col = 1 + labels, xlab = paste(c("axis_"),
axis_ind[1], sep = ""), ylab = paste(c("axis_"),
axis_ind[2], sep = ""))
points(pred_centers[[T]][, axis_ind[1]], pred_centers[[T]][,
axis_ind[2]], pch = 17, col = 1)
lgd = c("pred_centers ", sapply(1:length(unique(labels)),
function(x) paste(c("cluster"), as.character(x),
sep = " "), simplify = "vector"))
legend("topright", legend = lgd, col = seq(length(unique(labels)) +
1), pch = c(17, rep(1, length(unique(labels)))),
cex = 0.7, xjust = 1)
}
return(list(centers = pred_centers[[T]], nb_of_clusters = nb_of_clusters[T],
labels = labels))
}
else {
lambda_1 = coeff * (d + 2) * (T - 1)^(-0.5)/R
tau_proposal = (K_max * T * d)^(-0.5)
Nclusters[1] = 1
parameter_means_proposal[[1]] = kmeans(matrix(mydata,
ncol = d), Nclusters[1], nstart = 2, iter.max = 10)$centers
Niter_centers[[1]] = matrix(rmt(1, mean = as.vector(parameter_means_proposal[[1]]),
S = diag(tau_proposal, d), df = 3), nrow = 1)
while (sqrt(sum((Niter_centers[[1]]^2))) > multiplier_R *
R) {
Niter_centers[[1]] = matrix(rmt(1, mean = as.vector(parameter_means_proposal[[1]]),
S = diag(tau_proposal, d), df = 3), nrow = 1)
}
sum_loss[1] = cumulative_loss(Niter_centers[[1]],
mydata)
for (n in 2:N_iterations) {
transition_prob = c(transition_probability(Nclusters[n -
1], Nclusters[n - 1] - 1, K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1], K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1] + 1, K_max))
transition_prob = transition_prob/sum(transition_prob)
new_k = sample(c(Nclusters[n - 1] - 1, Nclusters[n -
1], Nclusters[n - 1] + 1), 1, prob = transition_prob)
if (new_k == Nclusters[n - 1]) {
m_t = parameter_means_proposal[[n - 1]]
}
else {
if (new_k >= T) {
new_k = T
m_t = matrix(mydata, ncol = d, byrow = F)
}
else {
m_t = kmeans(matrix(mydata, ncol = d), new_k,
nstart = 2, iter.max = 50)$centers
}
}
c_k_prime = t(apply(m_t, 1, function(x) rmt(1,
mean = x, S = diag(tau_proposal, d), df = 3)))
if (sum(sqrt(rowSums(c_k_prime^2)) < multiplier_R *
R) == new_k) {
log_numerator_prop = log(apply(matrix(1:Nclusters[n -
1], ncol = 1), 1, function(x) dmt(Niter_centers[[n -
1]][x, ], mean = parameter_means_proposal[[n -
1]][x, ], S = diag(tau_proposal, d), df = 3)))
log_denominator_prop = log(apply(matrix(1:new_k,
ncol = 1), 1, function(x) dmt(c_k_prime[x,
], mean = m_t[x, ], S = diag(tau_proposal,
d), df = 3)))
log_numerator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), new_k)
log_denominator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), Nclusters[n - 1])
ln_division = sum(c(log_numerator_prior, log_numerator_prop) -
c(log_denominator_prior, log_denominator_prop))
s_loss_prime = cumulative_loss(c_k_prime, mydata)
ln_accept_ratio = (-lambda_1 * (s_loss_prime -
sum_loss[n - 1])) + ln_division + log(transition_probability(new_k,
Nclusters[n - 1], K_max)) - log(transition_probability(Nclusters[n -
1], new_k, K_max))
}
else {
ln_accept_ratio = log(0)
}
bool = (log(runif(1)) < ln_accept_ratio)
if (bool) {
Niter_centers[[n]] = c_k_prime
parameter_means_proposal[[n]] = m_t
Nclusters[n] = new_k
sum_loss[n] = s_loss_prime
}
else {
Niter_centers[[n]] = Niter_centers[[n - 1]]
parameter_means_proposal[[n]] = parameter_means_proposal[[n -
1]]
Nclusters[n] = Nclusters[n - 1]
sum_loss[n] = sum_loss[n - 1]
}
}
pred_centers = Niter_centers[[N_iterations]]
nb_clusters = Nclusters[N_iterations]
labels = labels_function(pred_centers, mydata)
if (plot_ind) {
plot(mydata[, axis_ind[1]], mydata[, axis_ind[2]],
xlim = c(min(mydata[, axis_ind[1]]) - 5, max(mydata[,
axis_ind[1]]) + 5), ylim = c(min(mydata[,
axis_ind[2]]) - 5, max(mydata[, axis_ind[2]]) +
5), col = 1 + labels, xlab = paste(c("axis_"),
axis_ind[1], sep = ""), ylab = paste(c("axis_"),
axis_ind[2], sep = ""))
points(pred_centers[, axis_ind[1]], pred_centers[,
axis_ind[2]], col = 1, pch = 17)
lgd = c("pred_centers ", sapply(1:length(unique(labels)),
function(x) paste(c("cluster"), as.character(x),
sep = " "), simplify = "vector"))
legend("topright", legend = lgd, col = seq(length(unique(labels)) +
1), pch = c(17, rep(1, length(unique(labels)))),
cex = 0.7, xjust = 1)
}
return(list(predicted_centers = pred_centers, nb_of_clusters = nb_clusters,
labels = labels))
}
}
}
<bytecode: 0x244f6f8>
<environment: namespace:PACBO>
--- function search by body ---
Function PACBO in namespace PACBO has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(mydata) != "matrix") { : the condition has length > 1
Calls: PACBO
Execution halted
Flavor: r-devel-linux-x86_64-debian-clang
Version: 0.1.0
Check: examples
Result: ERROR
Running examples in ‘PACBO-Ex.R’ failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: PACBO
> ### Title: PACBO
> ### Aliases: PACBO
>
> ### ** Examples
>
> ## generating 4 clusters of 100 points in \strong{R}^{5}.
> set.seed(100)
> Nb <- 4
> d <- 5
> T <- 100
> proportion = rep(1/Nb, Nb)
> Mean_vectors <- matrix(runif(d*Nb,min=-10, max=10),nrow=Nb,ncol=d, byrow=TRUE)
> mydata <- matrix(replicate(T, rmnorm(1, mean= Mean_vectors[sample(1:Nb, 1, prob = proportion),],
+ varcov = diag(1,d))), nrow = T, byrow=T)
> R <- max(sqrt(rowSums(mydata^2)))
> ##run the algorithm.
> result <- PACBO(mydata, R, plot_ind = TRUE)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
PACBO
--- call from context ---
PACBO(mydata, R, plot_ind = TRUE)
--- call from argument ---
if (class(mydata) != "matrix") {
mydata = matrix(mydata, nrow = 1)
}
--- R stacktrace ---
where 1: PACBO(mydata, R, plot_ind = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE TRUE
--- function from context ---
function (mydata, R, coeff = 2, K_max = 50, scaling = FALSE,
var_ind = FALSE, N_iterations = 500, plot_ind = FALSE, axis_ind = c(1,
2))
{
if (class(mydata) != "matrix") {
mydata = matrix(mydata, nrow = 1)
}
if (R < max(sqrt(rowSums(mydata^2)))) {
print(c("R should be bigger than the maximum Euclidean distance of observations"))
}
else {
if (scaling) {
mydata = scale(mydata)
}
d = length(mydata[1, ])
T = length(mydata[, 1])
multiplier_R = 1.5
Nclusters = rep(1, N_iterations)
Niter_centers = list()
parameter_means_proposal = list()
sum_loss = rep(0, N_iterations)
if (var_ind) {
nb_of_clusters = rep(1, T)
pred_centers = list()
predicted_loss = rep(0, T)
proposal_loss = rep(0, T)
lambda_1 = rep(1, T)
lambda_1[2:T] = sapply(2:T, function(t) 2.5 * coeff *
(d + 2) * (t - 1)^(-0.5)/R)
lambda_2 = sapply(1:T, function(t) (d + 2) * (t)^(-0.5)/(multiplier_R *
R)^4)
c_1 = runiform_ball(5000, d, multiplier_R * R)
index_1 = which(apply((c_1 - t(replicate(5000, mydata[1,
])))^2, 1, sum) == instantaneous_loss(c_1, mydata[1,
]))
pred_centers[[1]] = c_1[index_1, ]
predicted_loss[1] = instantaneous_loss(pred_centers[[1]],
mydata[1, ])
pb <- txtProgressBar(0, 1, char = "=")
for (t in 2:T) {
tau_proposal = (K_max * t * d)^(-0.5)
Nclusters[1] = nb_of_clusters[t - 1]
parameter_means_proposal[[1]] = kmeans(matrix(mydata[1:(t -
1), ], ncol = d), Nclusters[1], nstart = 2,
iter.max = 10)$centers
Niter_centers[[1]] = t(apply(parameter_means_proposal[[1]],
1, function(x) rmt(1, mean = x, S = diag(tau_proposal,
d), df = 3)))
while (sum(sqrt(rowSums((Niter_centers[[1]]^2))) >
multiplier_R * R) > 0) {
Niter_centers[[1]] = t(apply(parameter_means_proposal[[1]],
1, function(x) rmt(1, mean = x, S = diag(tau_proposal,
d), df = 3)))
}
proposal_loss[1:(t - 1)] = apply(matrix(mydata[1:(t -
1), ], nrow = t - 1), 1, function(x) instantaneous_loss(Niter_centers[[1]],
x))
sum_loss[1] = sum(proposal_loss[1:(t - 1)]) +
0.5 * sum(lambda_2 * (proposal_loss - predicted_loss)^2)
for (n in 2:N_iterations) {
proposal_loss_temp = rep(0, T)
transition_prob = c(transition_probability(Nclusters[n -
1], Nclusters[n - 1] - 1, K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1], K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1] + 1, K_max))
transition_prob = transition_prob/sum(transition_prob)
new_k = sample(c(Nclusters[n - 1] - 1, Nclusters[n -
1], Nclusters[n - 1] + 1), 1, prob = transition_prob)
if (new_k == Nclusters[n - 1]) {
m_t = parameter_means_proposal[[n - 1]]
}
else {
if (new_k >= t - 1) {
new_k = t - 1
m_t = matrix(mydata[1:(t - 1), ], ncol = d,
byrow = F)
}
else {
m_t = kmeans(matrix(mydata[1:(t - 1), ],
ncol = d), new_k, nstart = 2, iter.max = 10)$centers
}
}
c_k_prime = t(apply(m_t, 1, function(x) rmt(1,
mean = x, S = diag(tau_proposal, d), df = 3)))
if (sum(sqrt(rowSums(c_k_prime^2)) < multiplier_R *
R) == new_k) {
log_numerator_prop = log(apply(matrix(1:Nclusters[n -
1], ncol = 1), 1, function(x) dmt(Niter_centers[[n -
1]][x, ], mean = parameter_means_proposal[[n -
1]][x, ], S = diag(tau_proposal, d), df = 3)))
log_denominator_prop = log(apply(matrix(1:new_k,
ncol = 1), 1, function(x) dmt(c_k_prime[x,
], mean = m_t[x, ], S = diag(tau_proposal,
d), df = 3)))
log_numerator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), new_k)
log_denominator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), Nclusters[n - 1])
ln_division = sum(c(log_numerator_prior,
log_numerator_prop) - c(log_denominator_prior,
log_denominator_prop))
proposal_loss_temp[1:(t - 1)] = apply(matrix(mydata[1:(t -
1), ], nrow = t - 1), 1, function(x) instantaneous_loss(c_k_prime,
x))
s_loss_prime = sum(proposal_loss_temp) +
0.5 * sum(lambda_2 * (proposal_loss_temp -
predicted_loss)^2)
ln_accept_ratio = (-lambda_1[t - 1] * (s_loss_prime -
sum_loss[n - 1])) + ln_division + log(transition_probability(new_k,
Nclusters[n - 1], K_max)) - log(transition_probability(Nclusters[n -
1], new_k, K_max))
}
else {
ln_accept_ratio = log(0)
}
bool = (log(runif(1)) < ln_accept_ratio)
if (bool) {
Niter_centers[[n]] = c_k_prime
parameter_means_proposal[[n]] = m_t
Nclusters[n] = new_k
sum_loss[n] = s_loss_prime
}
else {
Niter_centers[[n]] = Niter_centers[[n - 1]]
parameter_means_proposal[[n]] = parameter_means_proposal[[n -
1]]
Nclusters[n] = Nclusters[n - 1]
sum_loss[n] = sum_loss[n - 1]
}
}
nb_of_clusters[t] = Nclusters[N_iterations]
pred_centers[[t]] = Niter_centers[[N_iterations]]
predicted_loss[t] = instantaneous_loss(pred_centers[[t]],
mydata[t, ])
setTxtProgressBar(pb, t/T)
}
labels = labels_function(pred_centers[[T]], mydata)
if (plot_ind) {
plot(mydata[, axis_ind[1]], mydata[, axis_ind[2]],
xlim = c(min(mydata[, axis_ind[1]]) - 5, max(mydata[,
axis_ind[1]]) + 5), ylim = c(min(mydata[,
axis_ind[2]]) - 5, max(mydata[, axis_ind[2]]) +
5), col = 1 + labels, xlab = paste(c("axis_"),
axis_ind[1], sep = ""), ylab = paste(c("axis_"),
axis_ind[2], sep = ""))
points(pred_centers[[T]][, axis_ind[1]], pred_centers[[T]][,
axis_ind[2]], pch = 17, col = 1)
lgd = c("pred_centers ", sapply(1:length(unique(labels)),
function(x) paste(c("cluster"), as.character(x),
sep = " "), simplify = "vector"))
legend("topright", legend = lgd, col = seq(length(unique(labels)) +
1), pch = c(17, rep(1, length(unique(labels)))),
cex = 0.7, xjust = 1)
}
return(list(centers = pred_centers[[T]], nb_of_clusters = nb_of_clusters[T],
labels = labels))
}
else {
lambda_1 = coeff * (d + 2) * (T - 1)^(-0.5)/R
tau_proposal = (K_max * T * d)^(-0.5)
Nclusters[1] = 1
parameter_means_proposal[[1]] = kmeans(matrix(mydata,
ncol = d), Nclusters[1], nstart = 2, iter.max = 10)$centers
Niter_centers[[1]] = matrix(rmt(1, mean = as.vector(parameter_means_proposal[[1]]),
S = diag(tau_proposal, d), df = 3), nrow = 1)
while (sqrt(sum((Niter_centers[[1]]^2))) > multiplier_R *
R) {
Niter_centers[[1]] = matrix(rmt(1, mean = as.vector(parameter_means_proposal[[1]]),
S = diag(tau_proposal, d), df = 3), nrow = 1)
}
sum_loss[1] = cumulative_loss(Niter_centers[[1]],
mydata)
for (n in 2:N_iterations) {
transition_prob = c(transition_probability(Nclusters[n -
1], Nclusters[n - 1] - 1, K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1], K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1] + 1, K_max))
transition_prob = transition_prob/sum(transition_prob)
new_k = sample(c(Nclusters[n - 1] - 1, Nclusters[n -
1], Nclusters[n - 1] + 1), 1, prob = transition_prob)
if (new_k == Nclusters[n - 1]) {
m_t = parameter_means_proposal[[n - 1]]
}
else {
if (new_k >= T) {
new_k = T
m_t = matrix(mydata, ncol = d, byrow = F)
}
else {
m_t = kmeans(matrix(mydata, ncol = d), new_k,
nstart = 2, iter.max = 50)$centers
}
}
c_k_prime = t(apply(m_t, 1, function(x) rmt(1,
mean = x, S = diag(tau_proposal, d), df = 3)))
if (sum(sqrt(rowSums(c_k_prime^2)) < multiplier_R *
R) == new_k) {
log_numerator_prop = log(apply(matrix(1:Nclusters[n -
1], ncol = 1), 1, function(x) dmt(Niter_centers[[n -
1]][x, ], mean = parameter_means_proposal[[n -
1]][x, ], S = diag(tau_proposal, d), df = 3)))
log_denominator_prop = log(apply(matrix(1:new_k,
ncol = 1), 1, function(x) dmt(c_k_prime[x,
], mean = m_t[x, ], S = diag(tau_proposal,
d), df = 3)))
log_numerator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), new_k)
log_denominator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), Nclusters[n - 1])
ln_division = sum(c(log_numerator_prior, log_numerator_prop) -
c(log_denominator_prior, log_denominator_prop))
s_loss_prime = cumulative_loss(c_k_prime, mydata)
ln_accept_ratio = (-lambda_1 * (s_loss_prime -
sum_loss[n - 1])) + ln_division + log(transition_probability(new_k,
Nclusters[n - 1], K_max)) - log(transition_probability(Nclusters[n -
1], new_k, K_max))
}
else {
ln_accept_ratio = log(0)
}
bool = (log(runif(1)) < ln_accept_ratio)
if (bool) {
Niter_centers[[n]] = c_k_prime
parameter_means_proposal[[n]] = m_t
Nclusters[n] = new_k
sum_loss[n] = s_loss_prime
}
else {
Niter_centers[[n]] = Niter_centers[[n - 1]]
parameter_means_proposal[[n]] = parameter_means_proposal[[n -
1]]
Nclusters[n] = Nclusters[n - 1]
sum_loss[n] = sum_loss[n - 1]
}
}
pred_centers = Niter_centers[[N_iterations]]
nb_clusters = Nclusters[N_iterations]
labels = labels_function(pred_centers, mydata)
if (plot_ind) {
plot(mydata[, axis_ind[1]], mydata[, axis_ind[2]],
xlim = c(min(mydata[, axis_ind[1]]) - 5, max(mydata[,
axis_ind[1]]) + 5), ylim = c(min(mydata[,
axis_ind[2]]) - 5, max(mydata[, axis_ind[2]]) +
5), col = 1 + labels, xlab = paste(c("axis_"),
axis_ind[1], sep = ""), ylab = paste(c("axis_"),
axis_ind[2], sep = ""))
points(pred_centers[, axis_ind[1]], pred_centers[,
axis_ind[2]], col = 1, pch = 17)
lgd = c("pred_centers ", sapply(1:length(unique(labels)),
function(x) paste(c("cluster"), as.character(x),
sep = " "), simplify = "vector"))
legend("topright", legend = lgd, col = seq(length(unique(labels)) +
1), pch = c(17, rep(1, length(unique(labels)))),
cex = 0.7, xjust = 1)
}
return(list(predicted_centers = pred_centers, nb_of_clusters = nb_clusters,
labels = labels))
}
}
}
<bytecode: 0x564c3c38ad18>
<environment: namespace:PACBO>
--- function search by body ---
Function PACBO in namespace PACBO has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(mydata) != "matrix") { : the condition has length > 1
Calls: PACBO
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 0.1.0
Check: examples
Result: ERROR
Running examples in ‘PACBO-Ex.R’ failed
The error most likely occurred in:
> ### Name: PACBO
> ### Title: PACBO
> ### Aliases: PACBO
>
> ### ** Examples
>
> ## generating 4 clusters of 100 points in \strong{R}^{5}.
> set.seed(100)
> Nb <- 4
> d <- 5
> T <- 100
> proportion = rep(1/Nb, Nb)
> Mean_vectors <- matrix(runif(d*Nb,min=-10, max=10),nrow=Nb,ncol=d, byrow=TRUE)
> mydata <- matrix(replicate(T, rmnorm(1, mean= Mean_vectors[sample(1:Nb, 1, prob = proportion),],
+ varcov = diag(1,d))), nrow = T, byrow=T)
> R <- max(sqrt(rowSums(mydata^2)))
> ##run the algorithm.
> result <- PACBO(mydata, R, plot_ind = TRUE)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
PACBO
--- call from context ---
PACBO(mydata, R, plot_ind = TRUE)
--- call from argument ---
if (class(mydata) != "matrix") {
mydata = matrix(mydata, nrow = 1)
}
--- R stacktrace ---
where 1: PACBO(mydata, R, plot_ind = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE TRUE
--- function from context ---
function (mydata, R, coeff = 2, K_max = 50, scaling = FALSE,
var_ind = FALSE, N_iterations = 500, plot_ind = FALSE, axis_ind = c(1,
2))
{
if (class(mydata) != "matrix") {
mydata = matrix(mydata, nrow = 1)
}
if (R < max(sqrt(rowSums(mydata^2)))) {
print(c("R should be bigger than the maximum Euclidean distance of observations"))
}
else {
if (scaling) {
mydata = scale(mydata)
}
d = length(mydata[1, ])
T = length(mydata[, 1])
multiplier_R = 1.5
Nclusters = rep(1, N_iterations)
Niter_centers = list()
parameter_means_proposal = list()
sum_loss = rep(0, N_iterations)
if (var_ind) {
nb_of_clusters = rep(1, T)
pred_centers = list()
predicted_loss = rep(0, T)
proposal_loss = rep(0, T)
lambda_1 = rep(1, T)
lambda_1[2:T] = sapply(2:T, function(t) 2.5 * coeff *
(d + 2) * (t - 1)^(-0.5)/R)
lambda_2 = sapply(1:T, function(t) (d + 2) * (t)^(-0.5)/(multiplier_R *
R)^4)
c_1 = runiform_ball(5000, d, multiplier_R * R)
index_1 = which(apply((c_1 - t(replicate(5000, mydata[1,
])))^2, 1, sum) == instantaneous_loss(c_1, mydata[1,
]))
pred_centers[[1]] = c_1[index_1, ]
predicted_loss[1] = instantaneous_loss(pred_centers[[1]],
mydata[1, ])
pb <- txtProgressBar(0, 1, char = "=")
for (t in 2:T) {
tau_proposal = (K_max * t * d)^(-0.5)
Nclusters[1] = nb_of_clusters[t - 1]
parameter_means_proposal[[1]] = kmeans(matrix(mydata[1:(t -
1), ], ncol = d), Nclusters[1], nstart = 2,
iter.max = 10)$centers
Niter_centers[[1]] = t(apply(parameter_means_proposal[[1]],
1, function(x) rmt(1, mean = x, S = diag(tau_proposal,
d), df = 3)))
while (sum(sqrt(rowSums((Niter_centers[[1]]^2))) >
multiplier_R * R) > 0) {
Niter_centers[[1]] = t(apply(parameter_means_proposal[[1]],
1, function(x) rmt(1, mean = x, S = diag(tau_proposal,
d), df = 3)))
}
proposal_loss[1:(t - 1)] = apply(matrix(mydata[1:(t -
1), ], nrow = t - 1), 1, function(x) instantaneous_loss(Niter_centers[[1]],
x))
sum_loss[1] = sum(proposal_loss[1:(t - 1)]) +
0.5 * sum(lambda_2 * (proposal_loss - predicted_loss)^2)
for (n in 2:N_iterations) {
proposal_loss_temp = rep(0, T)
transition_prob = c(transition_probability(Nclusters[n -
1], Nclusters[n - 1] - 1, K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1], K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1] + 1, K_max))
transition_prob = transition_prob/sum(transition_prob)
new_k = sample(c(Nclusters[n - 1] - 1, Nclusters[n -
1], Nclusters[n - 1] + 1), 1, prob = transition_prob)
if (new_k == Nclusters[n - 1]) {
m_t = parameter_means_proposal[[n - 1]]
}
else {
if (new_k >= t - 1) {
new_k = t - 1
m_t = matrix(mydata[1:(t - 1), ], ncol = d,
byrow = F)
}
else {
m_t = kmeans(matrix(mydata[1:(t - 1), ],
ncol = d), new_k, nstart = 2, iter.max = 10)$centers
}
}
c_k_prime = t(apply(m_t, 1, function(x) rmt(1,
mean = x, S = diag(tau_proposal, d), df = 3)))
if (sum(sqrt(rowSums(c_k_prime^2)) < multiplier_R *
R) == new_k) {
log_numerator_prop = log(apply(matrix(1:Nclusters[n -
1], ncol = 1), 1, function(x) dmt(Niter_centers[[n -
1]][x, ], mean = parameter_means_proposal[[n -
1]][x, ], S = diag(tau_proposal, d), df = 3)))
log_denominator_prop = log(apply(matrix(1:new_k,
ncol = 1), 1, function(x) dmt(c_k_prime[x,
], mean = m_t[x, ], S = diag(tau_proposal,
d), df = 3)))
log_numerator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), new_k)
log_denominator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), Nclusters[n - 1])
ln_division = sum(c(log_numerator_prior,
log_numerator_prop) - c(log_denominator_prior,
log_denominator_prop))
proposal_loss_temp[1:(t - 1)] = apply(matrix(mydata[1:(t -
1), ], nrow = t - 1), 1, function(x) instantaneous_loss(c_k_prime,
x))
s_loss_prime = sum(proposal_loss_temp) +
0.5 * sum(lambda_2 * (proposal_loss_temp -
predicted_loss)^2)
ln_accept_ratio = (-lambda_1[t - 1] * (s_loss_prime -
sum_loss[n - 1])) + ln_division + log(transition_probability(new_k,
Nclusters[n - 1], K_max)) - log(transition_probability(Nclusters[n -
1], new_k, K_max))
}
else {
ln_accept_ratio = log(0)
}
bool = (log(runif(1)) < ln_accept_ratio)
if (bool) {
Niter_centers[[n]] = c_k_prime
parameter_means_proposal[[n]] = m_t
Nclusters[n] = new_k
sum_loss[n] = s_loss_prime
}
else {
Niter_centers[[n]] = Niter_centers[[n - 1]]
parameter_means_proposal[[n]] = parameter_means_proposal[[n -
1]]
Nclusters[n] = Nclusters[n - 1]
sum_loss[n] = sum_loss[n - 1]
}
}
nb_of_clusters[t] = Nclusters[N_iterations]
pred_centers[[t]] = Niter_centers[[N_iterations]]
predicted_loss[t] = instantaneous_loss(pred_centers[[t]],
mydata[t, ])
setTxtProgressBar(pb, t/T)
}
labels = labels_function(pred_centers[[T]], mydata)
if (plot_ind) {
plot(mydata[, axis_ind[1]], mydata[, axis_ind[2]],
xlim = c(min(mydata[, axis_ind[1]]) - 5, max(mydata[,
axis_ind[1]]) + 5), ylim = c(min(mydata[,
axis_ind[2]]) - 5, max(mydata[, axis_ind[2]]) +
5), col = 1 + labels, xlab = paste(c("axis_"),
axis_ind[1], sep = ""), ylab = paste(c("axis_"),
axis_ind[2], sep = ""))
points(pred_centers[[T]][, axis_ind[1]], pred_centers[[T]][,
axis_ind[2]], pch = 17, col = 1)
lgd = c("pred_centers ", sapply(1:length(unique(labels)),
function(x) paste(c("cluster"), as.character(x),
sep = " "), simplify = "vector"))
legend("topright", legend = lgd, col = seq(length(unique(labels)) +
1), pch = c(17, rep(1, length(unique(labels)))),
cex = 0.7, xjust = 1)
}
return(list(centers = pred_centers[[T]], nb_of_clusters = nb_of_clusters[T],
labels = labels))
}
else {
lambda_1 = coeff * (d + 2) * (T - 1)^(-0.5)/R
tau_proposal = (K_max * T * d)^(-0.5)
Nclusters[1] = 1
parameter_means_proposal[[1]] = kmeans(matrix(mydata,
ncol = d), Nclusters[1], nstart = 2, iter.max = 10)$centers
Niter_centers[[1]] = matrix(rmt(1, mean = as.vector(parameter_means_proposal[[1]]),
S = diag(tau_proposal, d), df = 3), nrow = 1)
while (sqrt(sum((Niter_centers[[1]]^2))) > multiplier_R *
R) {
Niter_centers[[1]] = matrix(rmt(1, mean = as.vector(parameter_means_proposal[[1]]),
S = diag(tau_proposal, d), df = 3), nrow = 1)
}
sum_loss[1] = cumulative_loss(Niter_centers[[1]],
mydata)
for (n in 2:N_iterations) {
transition_prob = c(transition_probability(Nclusters[n -
1], Nclusters[n - 1] - 1, K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1], K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1] + 1, K_max))
transition_prob = transition_prob/sum(transition_prob)
new_k = sample(c(Nclusters[n - 1] - 1, Nclusters[n -
1], Nclusters[n - 1] + 1), 1, prob = transition_prob)
if (new_k == Nclusters[n - 1]) {
m_t = parameter_means_proposal[[n - 1]]
}
else {
if (new_k >= T) {
new_k = T
m_t = matrix(mydata, ncol = d, byrow = F)
}
else {
m_t = kmeans(matrix(mydata, ncol = d), new_k,
nstart = 2, iter.max = 50)$centers
}
}
c_k_prime = t(apply(m_t, 1, function(x) rmt(1,
mean = x, S = diag(tau_proposal, d), df = 3)))
if (sum(sqrt(rowSums(c_k_prime^2)) < multiplier_R *
R) == new_k) {
log_numerator_prop = log(apply(matrix(1:Nclusters[n -
1], ncol = 1), 1, function(x) dmt(Niter_centers[[n -
1]][x, ], mean = parameter_means_proposal[[n -
1]][x, ], S = diag(tau_proposal, d), df = 3)))
log_denominator_prop = log(apply(matrix(1:new_k,
ncol = 1), 1, function(x) dmt(c_k_prime[x,
], mean = m_t[x, ], S = diag(tau_proposal,
d), df = 3)))
log_numerator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), new_k)
log_denominator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), Nclusters[n - 1])
ln_division = sum(c(log_numerator_prior, log_numerator_prop) -
c(log_denominator_prior, log_denominator_prop))
s_loss_prime = cumulative_loss(c_k_prime, mydata)
ln_accept_ratio = (-lambda_1 * (s_loss_prime -
sum_loss[n - 1])) + ln_division + log(transition_probability(new_k,
Nclusters[n - 1], K_max)) - log(transition_probability(Nclusters[n -
1], new_k, K_max))
}
else {
ln_accept_ratio = log(0)
}
bool = (log(runif(1)) < ln_accept_ratio)
if (bool) {
Niter_centers[[n]] = c_k_prime
parameter_means_proposal[[n]] = m_t
Nclusters[n] = new_k
sum_loss[n] = s_loss_prime
}
else {
Niter_centers[[n]] = Niter_centers[[n - 1]]
parameter_means_proposal[[n]] = parameter_means_proposal[[n -
1]]
Nclusters[n] = Nclusters[n - 1]
sum_loss[n] = sum_loss[n - 1]
}
}
pred_centers = Niter_centers[[N_iterations]]
nb_clusters = Nclusters[N_iterations]
labels = labels_function(pred_centers, mydata)
if (plot_ind) {
plot(mydata[, axis_ind[1]], mydata[, axis_ind[2]],
xlim = c(min(mydata[, axis_ind[1]]) - 5, max(mydata[,
axis_ind[1]]) + 5), ylim = c(min(mydata[,
axis_ind[2]]) - 5, max(mydata[, axis_ind[2]]) +
5), col = 1 + labels, xlab = paste(c("axis_"),
axis_ind[1], sep = ""), ylab = paste(c("axis_"),
axis_ind[2], sep = ""))
points(pred_centers[, axis_ind[1]], pred_centers[,
axis_ind[2]], col = 1, pch = 17)
lgd = c("pred_centers ", sapply(1:length(unique(labels)),
function(x) paste(c("cluster"), as.character(x),
sep = " "), simplify = "vector"))
legend("topright", legend = lgd, col = seq(length(unique(labels)) +
1), pch = c(17, rep(1, length(unique(labels)))),
cex = 0.7, xjust = 1)
}
return(list(predicted_centers = pred_centers, nb_of_clusters = nb_clusters,
labels = labels))
}
}
}
<bytecode: 0x336d950>
<environment: namespace:PACBO>
--- function search by body ---
Function PACBO in namespace PACBO has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(mydata) != "matrix") { : the condition has length > 1
Calls: PACBO
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 0.1.0
Check: examples
Result: ERROR
Running examples in ‘PACBO-Ex.R’ failed
The error most likely occurred in:
> ### Name: PACBO
> ### Title: PACBO
> ### Aliases: PACBO
>
> ### ** Examples
>
> ## generating 4 clusters of 100 points in \strong{R}^{5}.
> set.seed(100)
> Nb <- 4
> d <- 5
> T <- 100
> proportion = rep(1/Nb, Nb)
> Mean_vectors <- matrix(runif(d*Nb,min=-10, max=10),nrow=Nb,ncol=d, byrow=TRUE)
> mydata <- matrix(replicate(T, rmnorm(1, mean= Mean_vectors[sample(1:Nb, 1, prob = proportion),],
+ varcov = diag(1,d))), nrow = T, byrow=T)
> R <- max(sqrt(rowSums(mydata^2)))
> ##run the algorithm.
> result <- PACBO(mydata, R, plot_ind = TRUE)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
PACBO
--- call from context ---
PACBO(mydata, R, plot_ind = TRUE)
--- call from argument ---
if (class(mydata) != "matrix") {
mydata = matrix(mydata, nrow = 1)
}
--- R stacktrace ---
where 1: PACBO(mydata, R, plot_ind = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE TRUE
--- function from context ---
function (mydata, R, coeff = 2, K_max = 50, scaling = FALSE,
var_ind = FALSE, N_iterations = 500, plot_ind = FALSE, axis_ind = c(1,
2))
{
if (class(mydata) != "matrix") {
mydata = matrix(mydata, nrow = 1)
}
if (R < max(sqrt(rowSums(mydata^2)))) {
print(c("R should be bigger than the maximum Euclidean distance of observations"))
}
else {
if (scaling) {
mydata = scale(mydata)
}
d = length(mydata[1, ])
T = length(mydata[, 1])
multiplier_R = 1.5
Nclusters = rep(1, N_iterations)
Niter_centers = list()
parameter_means_proposal = list()
sum_loss = rep(0, N_iterations)
if (var_ind) {
nb_of_clusters = rep(1, T)
pred_centers = list()
predicted_loss = rep(0, T)
proposal_loss = rep(0, T)
lambda_1 = rep(1, T)
lambda_1[2:T] = sapply(2:T, function(t) 2.5 * coeff *
(d + 2) * (t - 1)^(-0.5)/R)
lambda_2 = sapply(1:T, function(t) (d + 2) * (t)^(-0.5)/(multiplier_R *
R)^4)
c_1 = runiform_ball(5000, d, multiplier_R * R)
index_1 = which(apply((c_1 - t(replicate(5000, mydata[1,
])))^2, 1, sum) == instantaneous_loss(c_1, mydata[1,
]))
pred_centers[[1]] = c_1[index_1, ]
predicted_loss[1] = instantaneous_loss(pred_centers[[1]],
mydata[1, ])
pb <- txtProgressBar(0, 1, char = "=")
for (t in 2:T) {
tau_proposal = (K_max * t * d)^(-0.5)
Nclusters[1] = nb_of_clusters[t - 1]
parameter_means_proposal[[1]] = kmeans(matrix(mydata[1:(t -
1), ], ncol = d), Nclusters[1], nstart = 2,
iter.max = 10)$centers
Niter_centers[[1]] = t(apply(parameter_means_proposal[[1]],
1, function(x) rmt(1, mean = x, S = diag(tau_proposal,
d), df = 3)))
while (sum(sqrt(rowSums((Niter_centers[[1]]^2))) >
multiplier_R * R) > 0) {
Niter_centers[[1]] = t(apply(parameter_means_proposal[[1]],
1, function(x) rmt(1, mean = x, S = diag(tau_proposal,
d), df = 3)))
}
proposal_loss[1:(t - 1)] = apply(matrix(mydata[1:(t -
1), ], nrow = t - 1), 1, function(x) instantaneous_loss(Niter_centers[[1]],
x))
sum_loss[1] = sum(proposal_loss[1:(t - 1)]) +
0.5 * sum(lambda_2 * (proposal_loss - predicted_loss)^2)
for (n in 2:N_iterations) {
proposal_loss_temp = rep(0, T)
transition_prob = c(transition_probability(Nclusters[n -
1], Nclusters[n - 1] - 1, K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1], K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1] + 1, K_max))
transition_prob = transition_prob/sum(transition_prob)
new_k = sample(c(Nclusters[n - 1] - 1, Nclusters[n -
1], Nclusters[n - 1] + 1), 1, prob = transition_prob)
if (new_k == Nclusters[n - 1]) {
m_t = parameter_means_proposal[[n - 1]]
}
else {
if (new_k >= t - 1) {
new_k = t - 1
m_t = matrix(mydata[1:(t - 1), ], ncol = d,
byrow = F)
}
else {
m_t = kmeans(matrix(mydata[1:(t - 1), ],
ncol = d), new_k, nstart = 2, iter.max = 10)$centers
}
}
c_k_prime = t(apply(m_t, 1, function(x) rmt(1,
mean = x, S = diag(tau_proposal, d), df = 3)))
if (sum(sqrt(rowSums(c_k_prime^2)) < multiplier_R *
R) == new_k) {
log_numerator_prop = log(apply(matrix(1:Nclusters[n -
1], ncol = 1), 1, function(x) dmt(Niter_centers[[n -
1]][x, ], mean = parameter_means_proposal[[n -
1]][x, ], S = diag(tau_proposal, d), df = 3)))
log_denominator_prop = log(apply(matrix(1:new_k,
ncol = 1), 1, function(x) dmt(c_k_prime[x,
], mean = m_t[x, ], S = diag(tau_proposal,
d), df = 3)))
log_numerator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), new_k)
log_denominator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), Nclusters[n - 1])
ln_division = sum(c(log_numerator_prior,
log_numerator_prop) - c(log_denominator_prior,
log_denominator_prop))
proposal_loss_temp[1:(t - 1)] = apply(matrix(mydata[1:(t -
1), ], nrow = t - 1), 1, function(x) instantaneous_loss(c_k_prime,
x))
s_loss_prime = sum(proposal_loss_temp) +
0.5 * sum(lambda_2 * (proposal_loss_temp -
predicted_loss)^2)
ln_accept_ratio = (-lambda_1[t - 1] * (s_loss_prime -
sum_loss[n - 1])) + ln_division + log(transition_probability(new_k,
Nclusters[n - 1], K_max)) - log(transition_probability(Nclusters[n -
1], new_k, K_max))
}
else {
ln_accept_ratio = log(0)
}
bool = (log(runif(1)) < ln_accept_ratio)
if (bool) {
Niter_centers[[n]] = c_k_prime
parameter_means_proposal[[n]] = m_t
Nclusters[n] = new_k
sum_loss[n] = s_loss_prime
}
else {
Niter_centers[[n]] = Niter_centers[[n - 1]]
parameter_means_proposal[[n]] = parameter_means_proposal[[n -
1]]
Nclusters[n] = Nclusters[n - 1]
sum_loss[n] = sum_loss[n - 1]
}
}
nb_of_clusters[t] = Nclusters[N_iterations]
pred_centers[[t]] = Niter_centers[[N_iterations]]
predicted_loss[t] = instantaneous_loss(pred_centers[[t]],
mydata[t, ])
setTxtProgressBar(pb, t/T)
}
labels = labels_function(pred_centers[[T]], mydata)
if (plot_ind) {
plot(mydata[, axis_ind[1]], mydata[, axis_ind[2]],
xlim = c(min(mydata[, axis_ind[1]]) - 5, max(mydata[,
axis_ind[1]]) + 5), ylim = c(min(mydata[,
axis_ind[2]]) - 5, max(mydata[, axis_ind[2]]) +
5), col = 1 + labels, xlab = paste(c("axis_"),
axis_ind[1], sep = ""), ylab = paste(c("axis_"),
axis_ind[2], sep = ""))
points(pred_centers[[T]][, axis_ind[1]], pred_centers[[T]][,
axis_ind[2]], pch = 17, col = 1)
lgd = c("pred_centers ", sapply(1:length(unique(labels)),
function(x) paste(c("cluster"), as.character(x),
sep = " "), simplify = "vector"))
legend("topright", legend = lgd, col = seq(length(unique(labels)) +
1), pch = c(17, rep(1, length(unique(labels)))),
cex = 0.7, xjust = 1)
}
return(list(centers = pred_centers[[T]], nb_of_clusters = nb_of_clusters[T],
labels = labels))
}
else {
lambda_1 = coeff * (d + 2) * (T - 1)^(-0.5)/R
tau_proposal = (K_max * T * d)^(-0.5)
Nclusters[1] = 1
parameter_means_proposal[[1]] = kmeans(matrix(mydata,
ncol = d), Nclusters[1], nstart = 2, iter.max = 10)$centers
Niter_centers[[1]] = matrix(rmt(1, mean = as.vector(parameter_means_proposal[[1]]),
S = diag(tau_proposal, d), df = 3), nrow = 1)
while (sqrt(sum((Niter_centers[[1]]^2))) > multiplier_R *
R) {
Niter_centers[[1]] = matrix(rmt(1, mean = as.vector(parameter_means_proposal[[1]]),
S = diag(tau_proposal, d), df = 3), nrow = 1)
}
sum_loss[1] = cumulative_loss(Niter_centers[[1]],
mydata)
for (n in 2:N_iterations) {
transition_prob = c(transition_probability(Nclusters[n -
1], Nclusters[n - 1] - 1, K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1], K_max), transition_probability(Nclusters[n -
1], Nclusters[n - 1] + 1, K_max))
transition_prob = transition_prob/sum(transition_prob)
new_k = sample(c(Nclusters[n - 1] - 1, Nclusters[n -
1], Nclusters[n - 1] + 1), 1, prob = transition_prob)
if (new_k == Nclusters[n - 1]) {
m_t = parameter_means_proposal[[n - 1]]
}
else {
if (new_k >= T) {
new_k = T
m_t = matrix(mydata, ncol = d, byrow = F)
}
else {
m_t = kmeans(matrix(mydata, ncol = d), new_k,
nstart = 2, iter.max = 50)$centers
}
}
c_k_prime = t(apply(m_t, 1, function(x) rmt(1,
mean = x, S = diag(tau_proposal, d), df = 3)))
if (sum(sqrt(rowSums(c_k_prime^2)) < multiplier_R *
R) == new_k) {
log_numerator_prop = log(apply(matrix(1:Nclusters[n -
1], ncol = 1), 1, function(x) dmt(Niter_centers[[n -
1]][x, ], mean = parameter_means_proposal[[n -
1]][x, ], S = diag(tau_proposal, d), df = 3)))
log_denominator_prop = log(apply(matrix(1:new_k,
ncol = 1), 1, function(x) dmt(c_k_prime[x,
], mean = m_t[x, ], S = diag(tau_proposal,
d), df = 3)))
log_numerator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), new_k)
log_denominator_prior = rep((log(gamma(d/2 +
1)) - (d/2) * log(pi) - d * log(multiplier_R *
R)), Nclusters[n - 1])
ln_division = sum(c(log_numerator_prior, log_numerator_prop) -
c(log_denominator_prior, log_denominator_prop))
s_loss_prime = cumulative_loss(c_k_prime, mydata)
ln_accept_ratio = (-lambda_1 * (s_loss_prime -
sum_loss[n - 1])) + ln_division + log(transition_probability(new_k,
Nclusters[n - 1], K_max)) - log(transition_probability(Nclusters[n -
1], new_k, K_max))
}
else {
ln_accept_ratio = log(0)
}
bool = (log(runif(1)) < ln_accept_ratio)
if (bool) {
Niter_centers[[n]] = c_k_prime
parameter_means_proposal[[n]] = m_t
Nclusters[n] = new_k
sum_loss[n] = s_loss_prime
}
else {
Niter_centers[[n]] = Niter_centers[[n - 1]]
parameter_means_proposal[[n]] = parameter_means_proposal[[n -
1]]
Nclusters[n] = Nclusters[n - 1]
sum_loss[n] = sum_loss[n - 1]
}
}
pred_centers = Niter_centers[[N_iterations]]
nb_clusters = Nclusters[N_iterations]
labels = labels_function(pred_centers, mydata)
if (plot_ind) {
plot(mydata[, axis_ind[1]], mydata[, axis_ind[2]],
xlim = c(min(mydata[, axis_ind[1]]) - 5, max(mydata[,
axis_ind[1]]) + 5), ylim = c(min(mydata[,
axis_ind[2]]) - 5, max(mydata[, axis_ind[2]]) +
5), col = 1 + labels, xlab = paste(c("axis_"),
axis_ind[1], sep = ""), ylab = paste(c("axis_"),
axis_ind[2], sep = ""))
points(pred_centers[, axis_ind[1]], pred_centers[,
axis_ind[2]], col = 1, pch = 17)
lgd = c("pred_centers ", sapply(1:length(unique(labels)),
function(x) paste(c("cluster"), as.character(x),
sep = " "), simplify = "vector"))
legend("topright", legend = lgd, col = seq(length(unique(labels)) +
1), pch = c(17, rep(1, length(unique(labels)))),
cex = 0.7, xjust = 1)
}
return(list(predicted_centers = pred_centers, nb_of_clusters = nb_clusters,
labels = labels))
}
}
}
<bytecode: 0x43d84c8>
<environment: namespace:PACBO>
--- function search by body ---
Function PACBO in namespace PACBO has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(mydata) != "matrix") { : the condition has length > 1
Calls: PACBO
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc