Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I have a vector of sequences with presence (1) and absence (0), from were I have calculated the first order Markov Process.

This is how the data looks like:


For thsi vector I have calculated:

  • Probability of Presence (P1) and Absence (P0)
  • Probability of having Presence followed by Presence (P11), Presence followed by Absence (P10), Absence followed by Presence (P01) and Absence followed by Absence (P00)

The transition probabilities were obtained by using a loop the checks for sequence of 2 values: to calculate P_00 I am using

   P_00 : dataset[j]==0 & dataset[j+1]==0 , etc.

These are the results:



The idea is to populate the NA with presence/absence sequences according to the Probability values obtained for the sequence. The problem is that for now I am not being able to find the best and more adequate process for this problem and I am a newbie with this type of problems. Even if there is already a package that can do this for me, I would prefer to understand how can I do it.

share|improve this question
This looks like a mark-recapture problem, where you're basically guessing $P_{00}$, either through ML or Bayesian models. –  Roman Luštrik Jan 16 '13 at 11:05
@A.R which package do already what do you want to do ? and if I undesrtand you try to find a (best)method to replace you NA according to probs? –  agstudy Jan 16 '13 at 12:07
@agstudy for example with the HiddenMarkov package I can predict the states (using the probabilities obtained with the "observed states") and compare the result with the observed states. However I don't want to loose the "observed states" and just want to "predict" the NA states. –  A.R Jan 16 '13 at 13:15
add comment

migrated from stackoverflow.com Jan 17 '13 at 17:05

This question came from our site for professional and enthusiast programmers.

1 Answer

up vote 2 down vote accepted

This is best and most transparently done in bayesian modelling. Bayesian inference works using MCMC (Markov Chain Monte Carlo) simulations, and that's exactly what you need. Example model code in bugs (not tested!):

dataset[1] ~ dbern(P1) # you must somehow solve the first element
for (i in 2:n) { # markov chain - define how each value depends on the previous one
    dataset[i] ~ dbern(p.presence[dataset[i - 1]])

dbern(p) stands for Bernoulli distribution with probability p. The array p.presence can be defined as (using P* variables as you defined them):

p.presence[0] = P01 / P0   # probability of presence, given absence in previous step
p.presence[1] = P11 / P1   # probability of presence, given presence in previous step

This way you can pass p.presence as input data (as in your example), but you could also let bugs to estimate it from the data!! (this is much more common and reasonable). You can of course get posterior distributions and MCMC simulation samples for the missing values (NA) in the dataset, and compute various other statistics on it.

share|improve this answer
I agree that MCMC is the way to go for this. But don't you also need to worry about transitions from the NAs to what was observed? I don't think your method "cares" if the third element of A.R.'s data set is a 1 or a 0, which seems wrong. In short, I think you need to add dependence on dataset[i + 1], not just dataset[i - 1] –  David J. Harris Jan 18 '13 at 4:45
@David, no! This is why bayesian modelling is so great. The NAs will be replaced by MCMC simulation samples using the model. Look at the rules I specified in the model (for loop) - this is exactly the markov chain principle, and that's all you have to care about! Don't be bothered with some details like missing values :-). We are in Bayesian world! You don't have to understand how it works at all, but if you want, read something about MCMC sampling or gibbs sampling... –  Curious Jan 18 '13 at 9:38
(+1) Ah, right, you're using BUGS. Yeah, it should take care of that for you. I was thinking of a different approach that only looked at the NAs. Nice answer. –  David J. Harris Jan 18 '13 at 18:23
Can't I do something like this: model { for (i in 1:N) { dataset[i] ~ dbern(p.presence[dataset1[i]]) dataset1[i]<-dataset[i]+1 } } list(N=10,dataset=c(0,0,1,1,0,NA,NA,0,NA,NA),p.presence=c(0.60,0.40)) –  A.R Jan 24 '13 at 12:41
@A.R. 1) and have you tried p.presence[dataset[i-1]+1]? And you mistakenly use i instead of i-1! 2) I don't know, I was at course with Marc Kery... maybe his book Introduction to WinBUGS for Ecologists could be of interest for you, mostly if you are ecologist.. –  Curious Jan 24 '13 at 13:41
show 2 more comments

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.