As you pointed out, to sample individual entries of the transition matrix from normal distribution is problematic as row sums can be larger than one.

One way is to run markovchainFit with bootstrap method to generate multiple transition matrix, which can be used in simulation.

B <- markovchainFit(data=c(A[1,], NA, A[2,], NA, A[3,], NA, A[4,], NA, A[5,]),
                 name="weather",
                 method='bootstrap',
                 nboot=1000)

#Check the row sums are indeed 1s
rowSums(B$bootStrapSamples[[1]])
#cloudy   rain  sunny 
#     1      1      1 

So, to simulate with both (1) and (2) uncertainties, we may first sample a transition matrix from normal distribution with mean and SE provided by the model. Then followed by your method for random Markov chain sampling.

As you pointed out, to sample individual entries of the transition matrix independently from normal distribution is problematic as row sums can be larger than one.

To stochastically generate a transition matrix, we can use bootstrap method provided in markovchainFit.

B <- markovchainFit(data=c(A[1,], NA, A[2,], NA, A[3,], NA, A[4,], NA, A[5,]),
                 name="weather",
                 method='bootstrap',
                 nboot=1000) #number of models generated

#Check the row sums are indeed 1s
rowSums(B$bootStrapSamples[[1]])
#cloudy   rain  sunny 
#     1      1      1 

So, to simulate with both (1) and (2) uncertainties, we may first sample a transition matrix from bootstrapping normal distribution with mean and SE provided by the model. And followed by your method for random Markov chain sampling.

(2) The uncertainty due to small sample size could be viewed as binomial distribution with small sample size n. In your example, the expected distribution for "sunny" at step 10 with 5 observation is B(n=5, p=0.30408) / 5.

for(j in 1:10000){
LL<-list()
for(i in 1:5){
  mcWeather2 <- mcWeather
  mcWeather2@transitionMatrix <- mcWeather2@transitionMatrix + 
                                 rnorm(9, sd=B$standardError) # add parameter uncertainty
  mcWeather2@transitionMatrix[mcWeather2@transitionMatrix>1] = 1
  mcWeather2@transitionMatrix[mcWeather2@transitionMatrix<0] = 0

  LL[[i]]<-data.frame(cat=rmarkovchain(n = 10, object = mcWeather2, t0 = "sunny",include.t0 = TRUE),index=i,time=1:11)  
}

EDIT

As you pointed out, to sample individual entries of the transition matrix from normal distribution is problematic as row sums can be larger than one.

One way is to run markovchainFit with bootstrap method to generate multiple transition matrix, which can be used in simulation.

B <- markovchainFit(data=c(A[1,], NA, A[2,], NA, A[3,], NA, A[4,], NA, A[5,]),
                 name="weather",
                 method='bootstrap',
                 nboot=1000)

#Check the row sums are indeed 1s
rowSums(B$bootStrapSamples[[1]])
#cloudy   rain  sunny 
#     1      1      1 

(2) The uncertainty due to small sample size could be viewed as binomial distribution with small sample size n. In your example, the expected distribution for "sunny" at step 10 with 5 observation is B(n=5, p=0.30408) / 5.

for(j in 1:10000){
LL<-list()
for(i in 1:5){
  mcWeather2 <- mcWeather
  #mcWeather2@transitionMatrix <- mcWeather2@transitionMatrix + 
  #                               rnorm(9, sd=B$standardError) # add parameter uncertainty
  #mcWeather2@transitionMatrix[mcWeather2@transitionMatrix>1] = 1
  #mcWeather2@transitionMatrix[mcWeather2@transitionMatrix<0] = 0

  #EDIT: using bootstrap method for parameter uncertainty
  mcWeather2@transitionMatrix <- B$bootStrapSamples[[sample(length(B$bootStrapSamples), 1)]]
  
  LL[[i]]<-data.frame(cat=rmarkovchain(n = 10, object = mcWeather2, t0 = "sunny",include.t0 = TRUE),index=i,time=1:11)  
}