Translating a function from R into Matlab -


i trying reproduce function in r in matlab. not have experience r difficult me understand of code. here r code:

# model calculations # ================== # # function determine length of stay df. # # t = time since hospital admission in days # age = age @ 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) }) 

i having difficulty in understanding following lines of code , mean: 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,]

here matlab code in have tried reproduce r code:

function ft = losdf(age, stroketype, dest)  % function 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); 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 

when plot output, not giving me same values in matlab in r. cannot identify going wrong; know values of pt correct. think may have made error in setting initialstatevec or in defining ft(t)?

if can give me advice on have went wrong appreciated.

this line temp.mat <- matrix(c(1,0,1,0,1-p,p),3,2,byrow=true) creates matrix of 3 rows , 2 columns , putting elements give him filling rows rows.

let's p = 0.5 , temp.mat :

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

this line initial.state.vec <- rep(0,7) creates vector of 7 0 :

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

and line initial.state.vec[c(type,type+1)] <- temp.mat[type,], according value of type : 1, 2 or 3, puts row indexed type of matrix in specific indexes of vector calculated according value of type. means :

for type = 1, take 1st row of matrix , put in 1st , 2nd index of vector. type = 2, take 2nd row of matrix , put in index 2 , 3 , type = 3, take 3rd row , put in index 3 , 4.

example 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 

hope translation.


Comments

Popular posts from this blog

Command prompt result in label. Python 2.7 -

javascript - How do I use URL parameters to change link href on page? -

amazon web services - AWS Route53 Trying To Get Site To Resolve To www -