Here is my code.
forw<-function(x,Pi,delta)
{n <- length(x)
m <- 2
nu<- matrix(0, nrow = n, ncol = m)
#Initialisation
nu[1, ] <- delta[,x[1]]/sum(delta[,x[1]])
########
for (i in 2:n)
{
nu[i, ] <- (nu[i - 1, ] %*% Pi[(2*x[i-1]-1):(2*x[i-1]),(2*x[i]-1):(2*x[i])])/sum(nu[i - 1, ] %*% Pi[(2*x[i-1]-1):(2*x[i-1]),(2*x[i]-1):(2*x[i])])
}
return(nu);
}
#####################
#DATA INPUT
x<-c(2,1,2,1,2,2,1,2,2,1,2,2,1,2,2,1,2,1,2,1,2,2,1,2,2,1,2,2,1,2,2,1)
p=0.3
q=1-p
lambda1=1
lambda2=0
mju1=0
mju2=1
Pi=matrix(c(p*lambda1,p*(1-lambda1),p*(1-lambda1),1+p*lambda1-2*p,q*mju1,p-q*mju1, q*(1-mju1),1+q*mju1-p-q,p*lambda2,q-p*lambda2,p*(1-lambda2),1+p*lambda2-q-p,q*mju2,q*(1-mju2),q*(1-mju2),1+q*mju2-2*q),nrow=4, byrow=T, dimnames = list(c("11","21","12","22"), c("11","21","12","22")))
delta=matrix(rep(0.25,4), nrow=2, byrow=T)
####################################
alpha=forw(x,Pi,delta)
As you know, in R, if there are several max values, R chooses the first one. Since in my example delta, the uniform initial distribution, leads to all alphas being the same, the function output should be reflected as a sequence of 1’s independently of the observations.
However, as you'll see, the function output does not consist of 1's only.
The problem arises because at some point an error starts occurring although in practice alphas should be the same.
E.g., the forward probabilities look as follows:
print(alpha, digits=12)
[,1] [,2]
[1,] 0.5 0.5
[2,] 0.5 0.5
[3,] 0.5 0.5
[4,] 0.5 0.5
[5,] 0.5 0.5
[6,] 0.5 0.5
[7,] 0.5 0.5
[8,] 0.5 0.5
[9,] 0.5 0.5
[10,] 0.5 0.5
[11,] 0.5 0.5
[12,] 0.5 0.5
[13,] 0.5 0.5
[14,] 0.5 0.5
[15,] 0.5 0.5
[16,] 0.5 0.5
[17,] 0.5 0.5
[18,] 0.5 0.5
[19,] 0.5 0.5
[20,] 0.5 0.5
[21,] 0.5 0.5
[22,] 0.5 0.5
[23,] 0.5 0.5
[24,] 0.5 0.5
[25,] 0.5 0.5
[26,] 0.5 0.5
[27,] 0.5 0.5
[28,] 0.5 0.5
[29,] 0.5 0.5
The columns look identical to 12 digits, however, the following output contradicts that. (In theoretical calculations all alphas are identical).
alpha[,1]-alpha[,2]
[1] 0.000000e+00 0.000000e+00 0.000000e+00
[4] 0.000000e+00 0.000000e+00 0.000000e+00
[7] 0.000000e+00 0.000000e+00 0.000000e+00
[10] 0.000000e+00 0.000000e+00 0.000000e+00
[13] 0.000000e+00 0.000000e+00 0.000000e+00
[16] 0.000000e+00 0.000000e+00 0.000000e+00
[19] -1.110223e-16 -3.330669e-16 -3.885781e-16
[22] -4.996004e-16 -6.106227e-16 -7.216450e-16
[25] -8.881784e-16 -1.110223e-15 -1.332268e-15
[28] -1.498801e-15 -1.609823e-15
Question: how can these differences be eliminated? Since they appear for the reason that is unclear to me I am now worried about the extent I can trust R outputs at all.
