File: fittingFunctions.R

package info (click to toggle)
r-cran-markovchain 0.8.5-4-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 2,060 kB
  • sloc: cpp: 2,854; sh: 13; makefile: 2
file content (759 lines) | stat: -rw-r--r-- 25,363 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
#' Function to generate a sequence of states from homogeneous Markov chains.
#' 
#' Provided any \code{markovchain} object, it returns a sequence of 
#' states coming from the underlying stationary distribution. 
#' 
#' @param n Sample size
#' @param markovchain \code{markovchain} object
#' @param t0 The initial state
#' @param include.t0 Specify if the initial state shall be used
#' @param useRCpp Boolean. Should RCpp fast implementation being used? Default is yes.
#' 
#' @details A sequence of size n is sampled.
#' 
#' @return A Character Vector
#' 
#' @references A First Course in Probability (8th Edition), Sheldon Ross, Prentice Hall 2010
#' 
#' @author Giorgio Spedicato
#' 
#' @seealso \code{\link{markovchainFit}}
#' 
#' @examples 
#' # define the markovchain object
#' statesNames <- c("a", "b", "c")
#' mcB <- new("markovchain", states = statesNames, 
#'    transitionMatrix = matrix(c(0.2, 0.5, 0.3, 0, 0.2, 0.8, 0.1, 0.8, 0.1), 
#'    nrow = 3, byrow = TRUE, dimnames = list(statesNames, statesNames)))
#' 
#' # show the sequence
#' outs <- markovchainSequence(n = 100, markovchain = mcB, t0 = "a")
#' 
#' @export

markovchainSequence <-function (n, markovchain, t0 = sample(markovchain@states, 1),
                               include.t0 = FALSE, useRCpp = TRUE) {
  
  # check whether given initial state is possible state or not
  if (!(t0 %in% markovchain@states))
    stop("Error! Initial state not defined")
  
  # call to cpp implmentation of markovchainSequence
  if (useRCpp) {
    return(.markovchainSequenceRcpp(n, markovchain, t0, include.t0))
  }
  
  # R implementation of the function
  
  # create a sequence of size n initially not initialized
  chain <- rep(NA,n)
  
  # initial state
  state <- t0
  
  # populate the sequence
  for (i in 1:n) {
    # row probabilty corresponding to the current state
    rowProbs <- markovchain@transitionMatrix[state, ]
    
    # select the next state
    outstate <- sample(size = 1, x = markovchain@states, prob = rowProbs)
    
    # store the new state
    chain[i] <- outstate
    
    # update the current state
    state <- outstate
  }
  
  # output
  out <- chain
  
  # whether to include initial state or not
  if (include.t0) {
    out <- c(t0, out)
  }
  
  return(out)
}


##################
# random sampler #
##################

# check if the subsequent states are included in the previous ones
# check the validity of non homogeneous markovchain list
# object is a list of markovchain object
.checkSequence <- function(object) {
  # assume non homogeneous markovchain list is valid
  out <- TRUE
  
  # list of one transition matrix implies valid
  if (length(object) == 1) {
    return(out) 
  }
  
  # if number of transition matrices are more than one  
  for (i in 2:length(object)) {
    
    # select the states which are reachable in one step
    if(object[[i - 1]]@byrow) {
      reachable <- (colSums(object[[i - 1]]@transitionMatrix) != 0)
    } else {
      reachable <- (rowSums(object[[i - 1]]@transitionMatrix) != 0)
    }
    
    # possible states in the previous markovchain object
    statesNm1 <- states(object[[i - 1]])[reachable]
    
    # possible states in the current markovchain object
    statesN <- states(object[[i]])
    
    # common states 
    intersection <- intersect(statesNm1, statesN)
    
    # condition to check whether statesNm1 is a subset of statesN or not
    if (setequal(intersection, statesNm1) == FALSE) {
      out <- FALSE
      break
    }
    
  }
  
  return(out)
}

#' Function to generate a sequence of states from homogeneous or non-homogeneous Markov chains.
#' 
#' Provided any \code{markovchain} or \code{markovchainList} objects, it returns a sequence of 
#' states coming from the underlying stationary distribution. 
#' 
#' @param n Sample size
#' @param object Either a \code{markovchain} or a \code{markovchainList} object
#' @param what It specifies whether either a \code{data.frame} or a \code{matrix} 
#'        (each rows represent a simulation) or a \code{list} is returned.
#' @param useRCpp Boolean. Should RCpp fast implementation being used? Default is yes.
#' @param parallel Boolean. Should parallel implementation being used? Default is yes.
#' @param num.cores Number of Cores to be used
#' @param ... additional parameters passed to the internal sampler
#' 
#' @details When a homogeneous process is assumed (\code{markovchain} object) a sequence is 
#' sampled of size n. When a non - homogeneous process is assumed,
#' n samples are taken but the process is assumed to last from the begin to the end of the 
#' non-homogeneous markov process.
#' 
#' @return Character Vector, data.frame, list or matrix
#' 
#' @references A First Course in Probability (8th Edition), Sheldon Ross, Prentice Hall 2010
#' 
#' @author Giorgio Spedicato
#' 
#' @note Check the type of input
#' 
#' @seealso \code{\link{markovchainFit}}, \code{\link{markovchainSequence}}
#' 
#' @examples 
#' # define the markovchain object
#' statesNames <- c("a", "b", "c")
#' mcB <- new("markovchain", states = statesNames, 
#'    transitionMatrix = matrix(c(0.2, 0.5, 0.3, 0, 0.2, 0.8, 0.1, 0.8, 0.1), 
#'    nrow = 3, byrow = TRUE, dimnames = list(statesNames, statesNames)))
#' 
#' # show the sequence
#' outs <- rmarkovchain(n = 100, object = mcB, what = "list")
#' 
#' 
#' #define markovchainList object
#' statesNames <- c("a", "b", "c")
#' mcA <- new("markovchain", states = statesNames, transitionMatrix = 
#'    matrix(c(0.2, 0.5, 0.3, 0, 0.2, 0.8, 0.1, 0.8, 0.1), nrow = 3, 
#'    byrow = TRUE, dimnames = list(statesNames, statesNames)))
#' mcB <- new("markovchain", states = statesNames, transitionMatrix = 
#'    matrix(c(0.2, 0.5, 0.3, 0, 0.2, 0.8, 0.1, 0.8, 0.1), nrow = 3, 
#'    byrow = TRUE, dimnames = list(statesNames, statesNames)))
#' mcC <- new("markovchain", states = statesNames, transitionMatrix = 
#'    matrix(c(0.2, 0.5, 0.3, 0, 0.2, 0.8, 0.1, 0.8, 0.1), nrow = 3, 
#'    byrow = TRUE, dimnames = list(statesNames, statesNames)))
#' mclist <- new("markovchainList", markovchains = list(mcA, mcB, mcC)) 
#' 
#' # show the list of sequence
#' rmarkovchain(100, mclist, "list")
#'      
#' @export

rmarkovchain <- function(n, object, what = "data.frame", useRCpp = TRUE, parallel = FALSE, num.cores = NULL, ...) {
  
  # check the class of the object
  if (class(object) == "markovchain") {
    out <- markovchainSequence(n = n, markovchain = object, useRCpp = useRCpp, ...)
    return(out)
  }
    
  if (class(object) == "markovchainList")
  {
    #######################################################
    if(useRCpp && !parallel) {
      
      # if include.t0 is not passed as extra argument then set include.t0 as false
      include.t0 <- list(...)$include.t0
      include.t0 <- ifelse(is.null(include.t0), FALSE, include.t0)
      
      # check whether initial state is passed or not
      t0 <- list(...)$t0
      if (is.null(t0)) t0 <- character()
      
      # call fast cpp function
      dataList <- .markovchainListRcpp(n, object@markovchains, include.t0, t0)
      
      # format in which results to be returned
      if (what == "data.frame") {
        out <- data.frame(iteration = dataList[[1]], values = dataList[[2]])
      }
        
      else {
        # output in matrix format
        # each row is an independent sequence
        out <- matrix(data = dataList[[2]], nrow = n, byrow = TRUE)
        
        # output in list format
        if (what == "list") {
          # outlist <- list()
          # for (i in 1:nrow(out))
          #  outlist[[i]] <- out[i, ]
          # out <- outlist
          out <- as.list(data.frame(t(out), stringsAsFactors = FALSE))
          out <- unname(out)
        }
      } 
      return(out)
    }
    ##########################################################
    if(useRCpp && parallel) {
      
      # Calculate the number of cores
      # It's not good to use all cores
      no_cores <- max(1,parallel::detectCores() - 1)
      
      # number of cores specified should be less than or equal to maximum cores available
      if((! is.null(num.cores))  && num.cores <= no_cores + 1 && num.cores >= 1) {
        no_cores <- num.cores
      }
      
      RcppParallel::setThreadOptions(no_cores)
      
      # if include.t0 is not passed as extra argument then set include.t0 as false
      include.t0 <- list(...)$include.t0
      include.t0 <- ifelse(is.null(include.t0), FALSE, include.t0)
      
      # check whether initial state is passed or not
      t0 <- list(...)$t0
      if (is.null(t0)) t0 <- character()
      
      dataList <- .markovchainSequenceParallelRcpp(object, n, include.t0, t0)
      
      if(what == "list") return(dataList)
      
      # dimension of matrix to be returned
      nrow <- length(dataList)
      ncol <- length(dataList[[1]])   
      
      if(what == "matrix") {
        out <- matrix(unlist(dataList), nrow = nrow, ncol = ncol, byrow = TRUE)
        # for(i in 1:nrow) out[i, ] <- dataList[[i]]
        return(out)
      }
      
      iteration <- unlist(lapply(1:nrow, rep, times = ncol))
      values <- unlist(dataList)
      
      # if what id data frame
      # for(i in 1:nrow) {
        # iteration <- c(iteration, rep(i, ncol))
        # values <- append(values, dataList[[i]])
      # }
      
      return(data.frame(iteration = iteration, values = values))
    }
    
    ##########################################################
    if(!useRCpp && parallel) {
      # if include.t0 is not passed as extra argument then set include.t0 as false
      include.t0 <- list(...)$include.t0
      include.t0 <- ifelse(is.null(include.t0), FALSE, include.t0)
      
      # check whether initial state is passed or not
      t0 <- list(...)$t0
      if (is.null(t0)) t0 <- character()
      
      dataList <- .markovchainSequenceParallel(n, object, t0, num.cores, include.t0)
      
      if(what == "list") return(dataList)
      
      # dimension of matrix to be returned
      nrow <- length(dataList)
      ncol <- length(dataList[[1]])   
      
      if(what == "matrix") {
        out <- matrix(nrow = nrow, ncol = ncol)
        for(i in 1:nrow) out[i, ] <- dataList[[i]]
        return(out)
      }
      
      iteration <- numeric()
      values <- character()
      
      # if what id data frame
      for(i in 1:nrow) {
        iteration <- append(iteration, rep(i, ncol))
        values <- append(values, dataList[[i]])
      }
      
      return(data.frame(iteration = iteration, values = values))
      
    }
    ##########################################################
    
    # store list of markovchain object in object
    object <- object@markovchains
    
    # check the validity of markovchainList object
    verify <- .checkSequence(object = object)
    
    # show warning if sequence is invalid
    if (!verify) {
      warning("Warning: some states in the markovchain sequences are not contained in the following states!")
    }
    
    # helper vector
    iteration <- numeric()
    values <- character()
    
    # create one sequence in each iteration
    for (i in 1:n) {
      
      # the first iteration may include initial state
      sampledValues <- markovchainSequence(n = 1, markovchain = object[[1]], ...)
      outIter <- rep(i, length(sampledValues))
      
      # number of markovchain objects are more than one
      if (length(object) > 1) {
        for (j in 2:length(object)) {
          pos2take <- length(sampledValues)
          # select new state of the sequence from the old state
          # t0 refers to the old state
          newVals <-markovchainSequence(n = 1, markovchain = object[[j]], t0 = sampledValues[pos2take]) 
          
          # update in every iteration
          outIter <- c(outIter, i)
          sampledValues <- c(sampledValues, newVals)
        }
      }
      
      # populate the helper vectors
      iteration <- c(iteration, outIter)
      values <- c(values, sampledValues)
    }
    
    # defining the output
    if (what == "data.frame") {
      out <- data.frame(iteration = iteration, values = values)
    } else {
      # ouput in matrix format
      out <- matrix(data = values, nrow = n, byrow = TRUE)
      
      # store each row of the matrix in the list
      if (what == 'list') {
        outlist <- list()
        for (i in 1:nrow(out))
          outlist[[i]] <- out[i, ]
        out <- outlist
      }
    }
  }
  
  return(out)
}

######################################################################

# helper function to calculate one sequence
.markovchainSPHelper <- function(x, t0, mclist, include.t0) {
  # number of transition matrices
  n <- length(mclist@markovchains)
  
  # take care of initial state
  vin <- 0
  if(include.t0) vin <- 1
  
  # a character vector to store a single sequence
  seq <- character(length = n + vin)
  
  if(length(t0) == 0) {
    stateNames <- mclist@markovchains[[1]]@states 
    t0 <- sample(x = stateNames, size = 1,  prob = rep(1 / length(stateNames), length(stateNames)))
  }
  
  if(include.t0) seq[1] <- t0
  
  invisible(lapply(seq_len(n),
    function(i)
    {
      stateNames <<- mclist@markovchains[[i]]@states 
        byRow <- mclist@markovchains[[i]]@byrow

        # check whether transition matrix follows row-wise or column-wise fashion
        if(byRow) prob <- mclist@markovchains[[i]]@transitionMatrix[which(stateNames == t0), ]
        else prob <- mclist@markovchains[[i]]@transitionMatrix[, which(stateNames == t0)]

        # initial state for the next transition matrix
        t0 <<- sample(x = stateNames, size = 1,  prob = prob)

        # populate the sequence vector
        seq[i+vin] <<- t0
  }
  ))
  
  return(seq)
}

# Function to generate a list of sequence of states in parallel from non-homogeneous Markov chains.
# 
# Provided any markovchainList object, it returns a list of sequence of states coming 
# from the underlying stationary distribution. 
# 
# @param  n Sample size
# @param object markovchainList object
# @param t0 Initial state
# @param num.cores Number of cores
#   

.markovchainSequenceParallel <- function(n, object,
                                        t0 = character(),
                                        num.cores = NULL, include.t0 = FALSE) {
  # check for the validity of non-uniform markov chain
  verify <- .checkSequence(object@markovchains)
  if (!verify) {
    warning("Warning: some states in the markovchain sequences are not contained in the following states!")
  }
    
  # Calculate the number of cores
  # It's not good to use all cores
  no_cores <- max(1,parallel::detectCores() - 1)
  
  # number of cores specified should be less than or equal to maximum cores available
  if((! is.null(num.cores))  && num.cores <= no_cores + 1 && num.cores >= 1) {
    no_cores <- num.cores
  }
  
  # Initiate cluster
  cl <- parallel::makeCluster(no_cores)
  
  # export the variables to be used in the helper function
  # parallel::clusterExport(cl, "t0")
  
  # export the variables to be used in the helper function
  mclist <- object
  # parallel::clusterExport(cl, "mclist")
   
  # list of n sequence
  listSeq <- tryCatch(parallel::parLapply(cl, 1:n, .markovchainSPHelper, t0, mclist, include.t0), 
                      error=function(e) e, warning=function(w) w)  
  
  # release the resources
  parallel::stopCluster(cl)
  
  return(listSeq)
}

######################################################################

# function to fit a DTMC with Laplacian Smoother
.mcFitLaplacianSmooth <- function(stringchar, byrow, laplacian = 0.01) {
  
  # every element of the matrix store the number of times jth state appears just
  # after the ith state
  origNum <- createSequenceMatrix(stringchar = stringchar, toRowProbs = FALSE)
  
  # add laplacian  to the sequence matrix
  # why? to avoid the cases where sum of row is zero
  newNum <- origNum + laplacian
  
  # store sum of each row  in the vector
  newSumOfRow <- rowSums(newNum)
  
  # helper matrix to convert frequency matrix to transition matrix
  newDen <- matrix(rep(newSumOfRow, length(newSumOfRow)), byrow = FALSE, ncol = length(newSumOfRow))
  
  # transition matrix
  transMatr <- newNum / newDen
  
  # create a markovchain object
  outMc <- new("markovchain", transitionMatrix = transMatr, name = "Laplacian Smooth Fit")

  # transpose the transition matrix
  if (!byrow) {
    outMc@transitionMatrix <- t(outMc@transitionMatrix)
    outMc@byrow <- FALSE
  }
  
  # wrap markovchain object in a list
  out <- list(estimate = outMc)
  return(out)
}

# function that return a Markov Chain from a given matrix of observations
# .matr2Mc <- function(matrData, laplacian = 0) {
#   
#   # number of columns in the input matrix  
#   nCols <- ncol(matrData)
#   
#   # an empty character vector to store names of possible states
#   uniqueVals <- character()
#   
#   # populate uniqueVals with names of states 
#   for(i in 1:nCols) {
#     uniqueVals <- union(uniqueVals, unique(as.character(matrData[,i]))) 
#   }
#   
#   # possible states in lexicographical order
#   uniqueVals <- sort(uniqueVals)
#   
#   # create a contingency matrix which store the number of times 
#   # jth state appear just after the ith state
#   contingencyMatrix <- matrix(rep(0, length(uniqueVals)^2), ncol = length(uniqueVals))
#   
#   # set the names of rows and columns
#   rownames(contingencyMatrix) <- colnames(contingencyMatrix) <- uniqueVals
#   
#   # fill the contingency matrix
#   for (i in 1:nrow(matrData)) {
#     for (j in 2:nCols) {
#       # state in the ith row and (j-1)th column
#       stateBegin <- as.character(matrData[i, j-1])
#       
#       # index of beginning state 
#       whichRow <- which(uniqueVals == stateBegin)
#       
#       # state in the ith row and jth column
#       stateEnd <- as.character(matrData[i, j])
#       
#       # index of ending state 
#       whichCols <- which(uniqueVals == stateEnd)
#       
#       # update the contingency matrix
#       contingencyMatrix[whichRow, whichCols] <- contingencyMatrix[whichRow, whichCols] + 1
#     }
#   }
#   
#   # add laplacian correction if needed
#   contingencyMatrix <- contingencyMatrix + laplacian
#   
#   # take care of rows with all entries 0
#   sumOfRows <- rowSums(contingencyMatrix) 
#   for(i in 1:length(sumOfRows)) {
#     if(sumOfRows[i] == 0) {
#       contingencyMatrix[i, ] <- 1
#       sumOfRows[i] <- length(sumOfRows)
#     }
#   }
#   
#   # get a transition matrix and a DTMC
#   transitionMatrix <- contingencyMatrix / sumOfRows
#   
#   # markov chain object to be returned
#   outMc <- new("markovchain", transitionMatrix = transitionMatrix)
#  
#   return(outMc)
# }


#' @title markovchainListFit
#' 
#' @description  Given a data frame or a matrix (rows are observations, by cols 
#' the temporal sequence), it fits a non - homogeneous discrete time markov chain 
#' process (storing row). In particular a markovchainList of size = ncol - 1 is obtained
#' estimating transitions from the n samples given by consecutive column pairs.
#' 
#' @param data Either a matrix or a data.frame or a list object.
#' @param laplacian Laplacian correction (default 0).
#' @param byrow Indicates whether distinc stochastic processes trajectiories are shown in distinct rows.
#' @param name Optional name.
#' 
#' @details If \code{data} contains \code{NAs} then the transitions containing \code{NA} will be ignored.
#' @return A list containing two slots:
#' estimate (the estimate)
#' name
#' 
#' @examples
#' 
#' # using holson dataset
#' data(holson)
#' # fitting a single markovchain
#' singleMc <- markovchainFit(data = holson[,2:12])
#' # fitting a markovchainList
#' mclistFit <- markovchainListFit(data = holson[, 2:12], name = "holsonMcList")
#' @export
markovchainListFit <- function(data, byrow = TRUE, laplacian = 0, name) {
  
  # check the format of input data
  if (!any(is.list(data),is.data.frame(data),is.matrix(data))) {
    stop("Error: data must be either a matrix or a data.frame or a list")
  }
  
  freqMatrixes <- list() 
  
  if(is.list(data) == TRUE) {
    markovchains <- list()
    # list of frquency matrix
    freqMatrixes <- .mcListFitForList(data)
    
  } else{
    # if input is data frame convert it to matrix
    if(is.data.frame(data)) {
      data <- unname(as.matrix(data))
    }
    
    # make the entries row wise if it is not
    if(!byrow) {
      data <- t(data) 
    }
    
    # number of columns in the matrix
    nCols <- ncol(data)
    
    # fit by columns
    freqMatrixes <- lapply(seq_len(nCols-1), function(i){
      # (i-1)th transition matrix for transition from (i-1)th state to ith state
      matrData <- data[, c(i, i+1)]
      matrData[1, ] <- as.character(matrData[1, ])
      # checking particular data for NA values.
      validTransition <- any(apply(matrData, 1, function(x){ !any(is.na(x)) }))
      
      if(validTransition)
        createSequenceMatrix(matrData, toRowProbs = FALSE, sanitize = TRUE)
    
    })
    
    freqMatrixes <- freqMatrixes[ !sapply(freqMatrixes, is.null) ]
  }
  
  if(length(freqMatrixes) == 0) {
    return(list())
  }
  
  
  markovchains <- lapply(freqMatrixes, function(freqMatrix){
    # add laplacian correction
    freqMatrix <- freqMatrix + laplacian
    rSums <- rowSums(freqMatrix)
    # transition matrix
    tMatrix <- freqMatrix / rSums;
    
    estMc <- new("markovchain", transitionMatrix = tMatrix)
    estMc
  })
  
  # create markovchainList object
  outMcList <- new("markovchainList", markovchains = markovchains)
  
  # wrap the object in a list
  out <- list(estimate = outMcList)
  
  # set the name of markovchainList object as given in the argument
  if(!missing(name)) {
    out$estimate@name <- name 
  }
  
  return(out)
}  

#' A function to compute multinomial confidence intervals of DTMC
#' 
#' @description Return estimated transition matrix assuming a Multinomial Distribution
#' 
#' @param transitionMatrix An estimated transition matrix.
#' @param countsTransitionMatrix Empirical (conts) transition matrix, on which the \code{transitionMatrix} was performed.
#' @param confidencelevel confidence interval level.
#' 
#' @return Two matrices containing the confidence intervals.
#' 
#' @seealso \code{markovchainFit}
#' 
#' @references Constructing two-sided simultaneous confidence intervals 
#' for multinomial proportions for small counts in a large number of cells. 
#' Journal of Statistical Software 5(6) (2000)
#' 
#' @examples 
#' seq<-c("a", "b", "a", "a", "a", "a", "b", "a", "b", "a", "b", "a", "a", "b", "b", "b", "a")
#' mcfit<-markovchainFit(data=seq,byrow=TRUE)
#' seqmat<-createSequenceMatrix(seq)
#' multinomialConfidenceIntervals(mcfit$estimate@transitionMatrix, seqmat, 0.95)
#' @export
multinomialConfidenceIntervals<-function(transitionMatrix, countsTransitionMatrix, confidencelevel=0.95) {
  
  out<-.multinomialCIRcpp(transMat=transitionMatrix, seqMat=countsTransitionMatrix,confidencelevel=confidencelevel)
  return(out)

}


#' return a joint pdf of the number of visits to the various states of the DTMC
#' 
#' @description This function would return a joint pdf of the number of visits to
#' the various states of the DTMC during the first N steps.
#' 
#' @usage noofVisitsDist(markovchain,N,state)
#' 
#' @param markovchain a markovchain-class object
#' @param N no of steps
#' @param state the initial state
#' 
#' @details 
#' This function would return a joint pdf of the number of visits to
#' the various states of the DTMC during the first N steps.
#' 
#' @return a numeric vector depicting the above described probability density function.
#' 
#' @author Vandit Jain
#' 
#' @examples 
#' transMatr<-matrix(c(0.4,0.6,.3,.7),nrow=2,byrow=TRUE)
#' simpleMc<-new("markovchain", states=c("a","b"),
#'              transitionMatrix=transMatr, 
#'              name="simpleMc")   
#' noofVisitsDist(simpleMc,5,"a")
#' 
#' @export
noofVisitsDist <- function(markovchain,N = 5,state) {
  
  if(class(markovchain)!="markovchain")
    stop("please provide a valid markovchain-class object")
  
  if(N <= 0)
    stop("please enter positive number of steps")
  
  # the transition matrix
  Tmatrix <- markovchain@transitionMatrix
  
  # character vector of states of the markovchain
  stateNames <- states(markovchain)
  
  i<--1
  
  # initial state
  i <- which(stateNames == state)
  
  if(i==-1)
    stop("please provide a valid inital state")
  
  
  # call to Rcpp implementation of the function
  out <- .noofVisitsDistRCpp(Tmatrix,i,N)
  
  # adds state names names to the output vector
  names(out) <- stateNames
  out <- c(out)
  return(out)
  
}