Normality is not really a requirement for this simulation, monte carlo methods are independent of the distribution, but lets have a look of how the weight is distributed for the past 3 years.

# include ggplot2
library(ggplot2)

gghist = function (mydata, mycolname) {
  pl = ggplot(data = mydata)
  subvp = viewport(width=0.35, height=0.35, x=0.84, y=0.84)

  his = pl + 
        geom_histogram(aes_string(x=mycolname,y="..density.."),alpha=0.2) + 
        geom_density(aes_string(x=mycolname)) + 
        opts(title = names(mydata[mycolname]))

  qqp = pl + 
        geom_point(aes_string(sample=mycolname), stat="qq") + labs(x=NULL, y=NULL) + 
        opts(title = "QQ")

  print(his)
  print(qqp, vp = subvp)
}

gghist(mydata, "kg")

Its close to normal distribution and says my weight has been mostly(in average) 80.5kg.

predict.weight.in.days = function(days, inicialweight, deltavector) {
  weight = inicialweight
  for (i in 1:days) {
    weight = weight + sample(deltavector, 1, replace=TRUE)
  }
  return(weight)
}

# simulate it 10k times
mcWeightWeek = replicate(10000, predict.weight.in.days(7, lastweight, kgdelta))

summary(mcWeightWeek)


   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  70.80   74.35   75.05   75.05   75.75   78.80

Another good thing about monte carlo methods is that they give a distribution of the prediction, so its possible to get a feeling of how certain the average is; either very certain with a big central peak, or not that certain when the graph is flatter and all over the place:

gghist(data.frame(kg=mcWeightWeek), "kg")

days.to.weight = function(weight, inicialweight, deltavector) {
  target = inicialweight
  days = 0
  while (target > weight) {
    target = target + sample(deltavector, 1, replace=TRUE)
     days = days + 1
     if (days >= 1095) # if value too crazy just interrupt the loop
        break
  }
  return(days)
}

# simulate it 10k times
mcDays75 = replicate(10000, days.to.weight(75, lastweight, kgdelta))

summary(mcDays75)


   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.000   4.000   7.000   9.073  11.000  82.000

And the cumulative distribution:

# add dates to it, from today's date + #days
days75 = sort(Sys.Date() + mcDays75)

# get the ecdf values into a dataframe
days75.ecdf = summarize(data.frame(days=days75), days = unique(days), 
                        ecdf = ecdf(days)(unique(days)))

# date where its 85% sure i'll reach goal
prob85 = head(days75.ecdf[days75.ecdf$ecdf>0.85,],1)

# plot
ggplot(days75.ecdf, aes(days, ecdf)) + geom_step() +
       ylab("probability") + 
       geom_point(aes(x = prob85$days, y = prob85$ecdf)) +
       geom_text(aes(x = prob85$days, y = prob85$ecdf, 
                    label = paste("85% sure to be 75kg on",
                            format(prob85$days, "%a, %d %b %Y"))), 
                     hjust=-0.04)

Also note that, weight loss is faster at the beginning of a diet, it tends to slow down over time, so to keep the predictions valid we need to continue record the weight and re-run the predictions frequently.

But as you see the slow carb diet seems to work, even without exercise. Tim’s book is great, focusing on the smallest things possible for the bigger results(=efficiency).

Hard drive occupation prediction with R: part 1 and part 2, and thanks to Leandro Penz on the feedback.