Traduction d'une fonction de R dans Matlab

Question

J'essaie de reproduire une fonction dans R dans Matlab. Je n'ai pas beaucoup d'expérience avec R, donc il m'est difficile de comprendre une partie du code. Voici le code R:

# Model Calculations 
# ================== 
# 
# Function to determine length of stay df. 
# 
# t = time since hospital admission in days 
# age = age at admision in years 
# type = type of stroke (1. Haemorhagic, 2. Cerebral Infarction, 3. TIA) 
# 
los.df <- Vectorize(function(t, age, type,dest){ 
    par <- c(6.63570, -0.03652, -3.06931, 0.07153, -8.66118, 
    0.08801, 22.10156, 2.48820, 1.56162, 1.27849, 
    11.76860, 3.41989, 63.92514) 
    alpha1 <- par[1] 
    beta1 <- par[2] 
    alpha2 <- par[3] 
    beta2 <- par[4] 
    theta0 <- par[5] 
    theta1 <- par[6] 
    mu1 <- par[7] 
    mu2 <- par[8] 
    mu3 <- par[9] 
    mu4 <- 0 
    nu1 <- 0 
    nu2 <- 0 
    nu3 <- par[10] 
    nu4 <- 0 
    rho1 <- 0 
    rho2 <- par[11] 
    rho3 <- par[12] 
    rho4 <- par[13] 
    # 
    p <- exp(-exp(theta0 + theta1*age)) 
    temp.mat <- matrix(c(1,0,1,0,1-p,p),3,2,byrow=TRUE) 
    initial.state.vec <- rep(0,7) 
    initial.state.vec[c(type,type+1)] <- temp.mat[type,] 
    # 
    lambda1 <- exp(alpha1 + beta1*age) 
    lambda2 <- exp(alpha2 + beta2*age) 
    Q <- matrix(0,7,7) 
    Q[1,] <- c(-(lambda1+mu1+nu1+rho1), lambda1, 0, 0, mu1, nu1, rho1) 
    Q[2,] <- c(0, -(lambda2+mu2+nu2+rho2), lambda2, 0, mu2, nu2, rho2) 
    Q[3,] <- c(0, 0, -(mu3+nu3+rho3), 0, mu3, nu3, rho3) 
    Q[4,] <- c(0, 0, 0, -(mu4+nu4+rho4), mu4, nu4, rho4) 
    Pt <- expm(t/365*Q) 
    Ft <- sum(as.numeric(initial.state.vec %*% Pt)[dest:dest]) 
    return(Ft) 
}) 

J'éprouve des difficultés à comprendre les lignes suivantes de code et ce qu'ils signifient: temp.mat <- matrix(c(1,0,1,0,1-p,p),3,2,byrow=TRUE) initial.state.vec <- rep(0,7) initial.state.vec[c(type,type+1)] <- temp.mat[type,]

Voici mon code Matlab dans lequel j'ai essayé de reproduire le code R :

function Ft = losdf(age, strokeType, dest) 

% function to determine length of stay in hospitaal of stroke patients 
% t = time since admission (days); 
% age = age of patient; 
% strokeType = 1. Haemorhagic, 2. Cerebral Infarction, 3. TIA; 
% dest = 5.Death 6.Nursing Home 7. Usual Residence; 

alpha1 = 6.63570; 
beta1 = -0.03652; 
alpha2 = -3.06931; 
beta2 = 0.07153; 
theta0 = -8.66118; 
theta1 = 0.08801; 
mu1 = 22.10156; 
mu2 = 2.48820; 
mu3 = 1.56162; 
mu4 = 0; 
nu1 = 0; 
nu2 = 0; 
nu3 = 1.27849; 
nu4 = 0; 
rho1 = 0; 
rho2 = 11.76860; 
rho3 = 3.41989; 
rho4 = 63.92514; 

Ft = zeros(365,1); 
for t = 1:1:365 
p = (exp(-exp(theta0 + (theta1.*age)))); 

if strokeType == 1 
    initialstatevec = [1 0 0 0 0 0 0]; 
elseif strokeType == 2 
    initialstatevec = [0 1 0 0 0 0 0]; 
else 
    initialstatevec = [0 0 (1-p) p 0 0 0]; 
end 

lambda1 = exp(alpha1 + (beta1.*age)); 
lambda2 = exp(alpha2 + (beta2.*age)); 

Q = [ -(lambda1+mu1+nu1+rho1) lambda1 0 0 mu1 nu1 rho1; 
0 -(lambda2+mu2+nu2+rho2) lambda2 0 mu2 nu2 rho2; 
0 0 -(mu3+nu3+rho3) 0 mu3 nu3 rho3; 
0 0 0 -(mu4+nu4+rho4) mu4 nu4 rho4; 
0 0 0 0 0 0 0; 
0 0 0 0 0 0 0; 
0 0 0 0 0 0 0]; 

Pt = expm(t./365.*Q); 
Pt = Pt(strokeType, dest); 
Ft(t) = sum(initialstatevec.*Pt); 

end 

Quand je tracer la sortie, il ne me donne pas les mêmes valeurs dans Matlab que dans R. Je ne peux pas identifier où je me trompe; Je sais que mes valeurs de Pt sont correctes. Je pense que j'ai peut-être fait une erreur en configurant initialstatevec ou en définissant Ft (t)?

Si quelqu'un peut me donner des conseils sur où je me suis trompé ce serait très apprécié.

Answer 1

Cette ligne temp.mat <- matrix(c(1,0,1,0,1-p,p),3,2,byrow=TRUE) crée une matrice de 3 lignes et 2 colonnes et met les éléments que vous lui donnez en remplissant les lignes par des lignes.

Disons que p = 0.5 et temp.mat est comme ceci:

 [,1] [,2] 
[1,] 1.0 0.0 
[2,] 1.0 0.0 
[3,] 0.5 0.5 

Cette ligne initial.state.vec <- rep(0,7) crée un vecteur de sept 0:

> initial.state.vec 
[1] 0 0 0 0 0 0 0 

Et cette ligne initial.state.vec[c(type,type+1)] <- temp.mat[type,], selon la valeur de type: 1 , 2 ou 3, met la ligne indexée par type de votre matrice dans des index spécifiques de vecteur qui sont également calculés ng à la valeur de type. Cela signifie que:

Pour type = 1, vous prenez la 1ère ligne de votre matrice et vous mettez dans le 1er et 2ème index de votre vecteur. Pour type = 2, vous prenez la 2ème ligne de votre matrice et vous le mettez dans l'index 2 et 3 et pour type = 3, vous prenez la 3ème ligne et que vous mettez dans l'index 3 et 4.

Exemple avec type = 3:

temp.mat[type,] 
[1] 0.5 0.5 

initial.state.vec 
[1] 0.0 0.0 0.5 0.5 0.0 0.0 0.0 

Espérons que cela aidera votre traduction.