Eine Funktion von R in Matlab übersetzen

Question

Ich versuche, eine Funktion in R in Matlab zu reproduzieren. Ich habe nicht viel Erfahrung mit R, deshalb ist es schwierig für mich, etwas vom Code zu verstehen. Hier ist der R-Code:

# 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) 
}) 

Ich habe Schwierigkeiten beim Verständnis der folgenden Codezeilen und was sie bedeuten: 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,]

Hier ist mein Matlab-Code, in dem ich habe versucht, den R-Code zu reproduzieren :

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 

Wenn ich die Ausgabe plotten, es gibt mir nicht die gleichen Werte in Matlab, wie es in R. ist, kann ich nicht erkennen, wo ich falsch werde; Ich weiß, dass meine Pt-Werte korrekt sind. Ich glaube, ich habe beim Einrichten von initialstatevec oder beim Definieren von Ft (t) einen Fehler gemacht.

Wenn jemand mir Rat geben kann, wo ich falsch gelaufen bin, würde dies sehr geschätzt werden.

Answer 1

Diese Zeile temp.mat <- matrix(c(1,0,1,0,1-p,p),3,2,byrow=TRUE) erstellt eine Matrix aus 3 Zeilen und 2 Spalten und setzt die Elemente, die Sie ihm geben, indem Sie Zeilen für Zeilen füllen.

Sagen wir p = 0.5 und temp.mat ist wie folgt:

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

Diese Linie initial.state.vec <- rep(0,7) einen Vektor von sieben 0 erstellt:

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

Und diese Linie initial.state.vec[c(type,type+1)] <- temp.mat[type,], entsprechend dem Wert von type: 1 , 2 oder 3, setzt die durch type indizierte Zeile Ihrer Matrix in spezifische Vektorindizes, die ebenfalls entsprechend berechnet werden ng auf den Wert type. Es bedeutet, dass:

Für type = 1, nehmen Sie die erste Zeile Ihrer Matrix und Sie in den ersten und zweiten Index Ihres Vektors. Für type = 2, nehmen Sie die zweite Zeile der Matrix und setzen Sie sie in den Index 2 und 3 und für type = 3, nehmen Sie die 3. Reihe und Sie setzen in den Index 3 und 4

Beispiel mit 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 

Ich hoffe, es wird für Ihre Übersetzung helfen.