Preciso encontrar o gradiente de uma função, para isso, estou usando a linha final do código abaixo, onde se lê: gr2 <- ifelse(y>0,df1bi+df2bi+df3bi+df4bi,dFTbi)
, note que ao compilar a primeira entrada deste código (df1bi+df2bi+df3bi+df4bi
), tem-se uma matriz de 6 linhas e 2 colunas, mesma coisa ao compilar a segunda entrada (dFTbi
), porém quando compilo a função completa (gr2
) recebo a primeira coluna da matriz (df1bi+df2bi+df3bi+df4bi
) como um vetor. Preciso receber de volta a matriz df1bi+df2bi+df3bi+df4bi
com as entradas que possui NaN
, correspondentes a y==0
substituidos pelos valores correspondentes na matriz dFTbi
. Alguém possui alguma sugestão? Segue o código completo!
n <- 6
beta0=0.5
beta1=1
beta2=2
gamma0=1
gamma1=2
gamma2=3
phi1 <- 3
phi2 <- 1
rho <- 0.2
X <- t(cbind(rbind(runif(n/3,0,1),rep(0,n/3),rep(0,n/3)),rbind(rep(0,n/3),runif(n/3,0,1),rep(0,n/3)),rbind(rep(0,n/3),rep(0,n/3),runif(n/3,0,1))))
mu1 <- exp(X%*%rbind(beta0,beta1,beta2))
W <- t(cbind(rbind(runif(n/3,0,1),rep(0,n/3),rep(0,n/3)),rbind(rep(0,n/3),runif(n/3,0,1),rep(0,n/3)),rbind(rep(0,n/3),rep(0,n/3),runif(n/3,0,1))))
mu2 <- exp(W%*%rbind(gamma0,gamma1,gamma2))
set.seed(333)
y <- sample(seq(0,10,length.out = n),n,replace = TRUE)
dens.normal <- function(x){
(1/sqrt(2*pi))*exp((-x^2)/2)
}
term0 <- ifelse(y>0,((y*(phi1+1)/(phi1*mu1))^(1/2)-((phi1*mu1)/(y*(phi1+1)))^(1/2)),0)
term1 <- exp((-phi1/4)*(term0)^2)
term2 <- ifelse(y>0,(((phi1+1)/(phi1*mu1*y))^(1/2)+((phi1*mu1)/((phi1+1)*(y^3)))^(1/2)),0)
term3 <- (1/(2*sqrt(2*pi)))*((phi1/2)^(1/2))
term4 <- (((phi2+1))/(mu2*(1-rho^2)))^(1/2)
term5 <- (phi2*mu2-phi2-1)/(sqrt(2)*(phi2+1))
term6 <- rho*((phi1)/(2*(1-rho^2)))^(1/2)
integrand <- ifelse(y>0,term4*term5+term6*term0,0)
term7 <- ifelse(y>0,pnorm(integrand),0)
term8 <- ((phi2/2)^(1/2))*(((phi2+1)/(phi2*mu2))^(1/2)-((phi2*mu2)/(phi2+1))^(1/2))
term9 <- ifelse(y>0,((y*(phi1+1))/(phi1*(mu1)))^(1/2)+((phi1*mu1)/(y*(phi1+1)))^(1/2),0)#Derivada de term0 em relação a beta0
term10 <- ifelse(y>0,((mu1*phi1)/((y^3)*(phi1+1)))^(1/2)-((phi1+1)/(y*phi1*mu1))^(1/2),0)#Derivada de term0 em relação a beta0
FT2 <- pnorm(term8)
p <- ncol(X)-1
j <- 2:ncol(X)
dFTb0 <- matrix(0,ncol = 1,nrow = n)
dFTbi <- matrix(1,ncol = p,nrow = n)
#f1=log(term1), f2=log(term2), f3=log(term3), f4=log(term7), f5=log(term8)
df1b0 <- ifelse(y>0,(phi1/4)*term0*term9,0)
df2b0 <- ifelse(y>0,(1/2)*(term10/term2),0)
df1bi <- as.vector(df1b0)*matrix(X[,j],ncol = p,nrow = n)
df2bi <- as.vector(df2b0)*matrix(X[,j],ncol = p,nrow = n)
df3bi <- matrix(0,ncol = p,nrow = n)
df4bi <- as.vector(dens.normal(term4*term5+term6*term0)*(term9/term7))*matrix(X[,j],ncol = p,nrow = n) #derivada de log(T2) em relação a beta[i], i>0
gr2 <- ifelse(y>0,df1bi+df2bi+df3bi+df4bi,dFTbi)