multiplicative algorithm for two factors model.
library(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutfun.agorithm = function(delta,x1.power, x2.power, d.or.F = 'd',criterion = 'D', func = '1', num.simul = 3000000){
#----FUNTION CREATE VECTOR V----#
k1 = sapply(x1.power, function(t) x1^t)
k2 = sapply(x2.power, function(t) x2^t)
V = t(k1*k2)
#----use loop to run algorithm----#
repeat{
# use D-criterion OR A OR c #
M.inv = switch(criterion,
D = solve(V%*%diag(p0)%*%t(V)),
A = solve(V%*%diag(p0)%*%t(V))%*%solve(V%*%diag(p0)%*%t(V)))
dj = diag(t(V)%*%M.inv%*%V)
z = switch(d.or.F,
d = dj,
`F` = dj - sum(p0*dj) )
fun.pj = switch(func,
`1` = exp(delta*z),
`2` = (exp(delta*z))/(1+(exp(delta*z))),
`3` = pnorm(delta*z),
`4` = log(exp(1) + delta*z),
`5` = 1.0001 - exp(-delta*z),
`6` = (1+sign(z)*z)^(sign(z)*delta))
p1 = (p0*fun.pj)/sum(p0*fun.pj)
p0=p1;
if(i >= num.simul) {warning(paste('it reaches limited simulated number \n Max.F =',
round(max(dj - sum(p1*dj)),6),'& delta =',delta,' & criterion =',criterion ))
r = data.frame(design.points = 'NoExist',
weight = 'NoExist',
condition.F = '>>zero',
simulation = 'TooMuch'
)
break
}
if(is.nan(max(dj - sum(p1*dj)))) {stop(paste('condition max of F is Infinite \n & delta =',delta,
'& function', func, '& criterion',criterion))}
if(max(dj - sum(p1*dj)) <= 10^-10){
r = data.frame( X1 =x1,X2=x2,
weight = p1,
dj = dj,
condition.F = (dj - sum(p1*dj)),
simulation = i)
#r = r[abs(r$condition.F) <0.0001 ,]
r = round(r,6)
break
}
i = i+1
}#repeat loop
return(as.matrix(r))
}# function.algorithm
##############################################
##############################################
p1 = c(0,1,0,1,2,0)
p2 = c(0,0,1,1,0,2)
x <- y <- seq(-1,1,by = 0.1)
t = expand.grid(x,y)
x1 = t[,2]; x2 = t[,1]
p0 = rep(1/length(x1),length(x1)); i = 0
#----------------------------------------#
d0 = fun.agorithm(0.1,p1,p2,'d','D',5);
d0[abs(d0[,5])<0.0001,]## X1 X2 weight dj condition.F simulation
## [1,] -1 -1 0.145791 6 0 3293
## [2,] -1 0 0.080161 6 0 3293
## [3,] -1 1 0.145791 6 0 3293
## [4,] 0 -1 0.080161 6 0 3293
## [5,] 0 0 0.096193 6 0 3293
## [6,] 0 1 0.080161 6 0 3293
## [7,] 1 -1 0.145791 6 0 3293
## [8,] 1 0 0.080161 6 0 3293
## [9,] 1 1 0.145791 6 0 3293#------------------------------------#
m = matrix(d0[,4], nrow = 21)
d = data.frame(x,y)
m2 = matrix(d0[,3], nrow = 21)
plot_ly(x = d$x,y = d$y, z = m)%>% add_surfaceplot_ly(x = d$x,y = d$y,z = m2)%>% add_surfacesecondly, sequential algorithm for two factors models
library(plotly)
var.func = function(power1, power2, x){
k1 = sapply(power1, function(t) x[,1]^t)
k2 = sapply(power2, function(t) x[,2]^t)
V = t(k1*k2)
p0 = rep(1/nrow(x), nrow(x))
M.inv = solve(V%*%diag(p0)%*%t(V))
dj = diag(t(V)%*%M.inv%*%V)
return(dj)
}
p1 = c(0,1,0,1,2,0)
p2 = c(0,0,1,1,0,2)
x <- y <- seq(-1,1,by = 0.1)
k = expand.grid(x,y)
func.2.seq = function(simulation,t){
for(i in 1:simulation){
vars = var.func(p1,p2,t)
added = t[which(vars - max(vars) ==0),]
t = rbind(t,added)
}
r = table(t)/ sum(table(t))
return(r)
}
r4 = round(func.2.seq(50,k),6)
index = which(r4 > 0.001, arr.ind = TRUE)
cbind(x[index[,2]],x[index[,1]], r4[r4>0.001])## [,1] [,2] [,3]
## [1,] -1 -1 0.126233
## [2,] -1 0 0.126233
## [3,] -1 1 0.126233
## [4,] 0 -1 0.063116
## [5,] 0 0 0.063116
## [6,] 0 1 0.063116
## [7,] 1 -1 0.126233
## [8,] 1 0 0.126233
## [9,] 1 1 0.126233d = data.frame(x1 = x,x2=y)
plot_ly(d,x = d$x1,y = d$x2, z = r4) %>% add_surface()Thirdly, multilicative algorithm for one factor model
fun.agorithm = function(delta,model.power =2, d.or.F = 'd',criterion = 'D', func = '1', num.simul = 3000000){
#----FUNTION CREATE VECTOR V----#
vec.V = function(j){
v = c()
C = matrix(c(numeric(j),1),ncol = 1)
while(j>=0){
v = rbind(x^(j),v)
j = j-1
}
return(v)
}
#--------------------------------#
V = vec.V(model.power)
#----use loop to run algorithm----#
repeat{
# use D-criterion OR A OR c #
M.inv = switch(criterion,
D = solve(V%*%diag(p0)%*%t(V), tol = 1e-18),
A = solve(V%*%diag(p0)%*%t(V),tol = 1e-18)%*%solve(V%*%diag(p0)%*%t(V),tol = 1e-18))
dj = diag(t(V)%*%M.inv%*%V)
z = switch(d.or.F,
d = dj,
`F` = dj - sum(p0*dj) )
fun.pj = switch(func,
`1` = exp(delta*z),
`2` = (exp(delta*z))/(1+(exp(delta*z))),
`3` = pnorm(delta*z),
`4` = log(exp(1) + delta*z),
`5` = 1.0001 - exp(-delta*z),
`6` = (1+sign(z)*z)^(sign(z)*delta))
p1 = (p0*fun.pj)/sum(p0*fun.pj)
p0=p1;
if(i >= num.simul) {warning(paste('it reaches limited simulated number \n Max.F =',
round(max(dj - sum(p1*dj)),6),'& delta =',delta,' & criterion =',criterion ))
r = data.frame(design.points = 'NoExist',
probabilities = 'NoExist',
condition.F = '>>zero',
simulation = 'TooMuch'
)
break
}
if(is.nan(max(dj - sum(p1*dj)))) {stop(paste('condition max of F is Infinite \n & delta =',delta,
'& function', func, '& criterion',criterion))}
if(max(dj - sum(p1*dj)) <= 10^-10){
r = data.frame( design.points =x, probabilities = p1,
dj = dj,
condition.F = (dj - sum(p1*dj)),
simulation = i)
# r = r[abs(r$condition.F) <0.001 ,]
r = round(r,6)
break
}
i = i+1
}#repeat loop
return(as.matrix(r))
}# function.algorithm
####################################
####################################
par(mfrow = c(1,3))
x = seq(-1,1, by = 0.1)
# define initial value of p0
p0 = rep(1/length(x),length(x)); i = 0
#---------#
d0 = fun.agorithm(0.1, model.power = 3)
d0[abs(d0[,4]) < 0.0001,][,1:3]## design.points probabilities dj
## [1,] -1.0 0.249529 4
## [2,] -0.5 0.112861 4
## [3,] -0.4 0.137611 4
## [4,] 0.4 0.137611 4
## [5,] 0.5 0.112861 4
## [6,] 1.0 0.249529 4ggplotly(ggplot(data.frame(x=d0[,1],y=d0[,2]),aes(x,y))+
geom_line(color = 'chocolate')+ geom_point(color = 'chocolate3'))## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`ggplotly(ggplot(data.frame(x=d0[,1],y=d0[,3]), aes(x,y))+
geom_line(color = 'chocolate')+ geom_point(color = 'chocolate3'))## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`sequential algorithm for one factor model
x = seq(-1,1,by = 0.1)
func.result = function(x, simulation, power, graphs = TRUE){
#--------------------------------------#
func.dj = function(power,design.points){
V = sapply(design.points, function(b) b^power)
p0 = rep(1/length(design.points), length(design.points))
M_inv = solve(V%*%diag(p0)%*%t(V))
dj = diag(t(V)%*%M_inv%*%V)
return(dj)
}
#--------------------------------------#
for(i in 1:simulation){
dj = round(func.dj(0:power,x),6)
added.points = x[which(!(dj - max(dj)))]
x = sort(c(x,added.points))
}
r = as.matrix(data.frame(round(table(x)/sum(table(x)),6)))
return(r)
}
r= func.result(x,25,3)
r[r[,2] > 0.001,]## x Freq
## [1,] "-1" "0.248182"
## [2,] "-0.5" "0.124091"
## [3,] "-0.4" "0.124091"
## [4,] "0.4" "0.124091"
## [5,] "0.5" "0.124091"
## [6,] "1" "0.248182"ggplotly(ggplot(data.frame(x=as.numeric(r[,1]),y=as.numeric(r[,2])),aes(x,y))+geom_line(color = 'chocolate')+ geom_point(color = 'chocolate3'))## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`