How to perform a Monte Carlo simulation?

How to do a Monte Carlo simulation for the attached problem?

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

Answer 1

[Image Analyst]

Here's a hint/start:

N = 100; % Whatever....
fullBottle = ones(1, N)
k = 50; % Whatever...
for dayNumber = 1 : k
  % some code involving randi() to extract a random pill
end

[Image Analyst]

If you still haven't solved it, post your code - what you have so far.

Once we've solved that, why don't we make it more fun by computing the life expectancy of the dog? Since the pills are meant to be given once per month (12 times per year) and you're dosing half a pill a day (182 times per year), the dog is clearly being overdosed and that may have some negative side effects. Assume a certain metabolism rate, and a certain lethality dose (perhaps you can find out the "LD50" amount for heartworm medicine) and find out the expected length of time before the dog dies. You might want to ask your professor how many dogs per year he goes through.

Answer 2

[Roger Stafford]

Assuming you are not required to use a Monte Carlo approach, I would think that the person assigning this problem would expect you to use Matlab to exactly calculate the requested expected number. It is something Matlab can do very handily if programmed correctly. The trouble with using a Monte Carlo method is its accuracy dependence on the sample size. Unless you use a truly enormous number of samples, you are bound to have an inaccurate estimate of the statistics, particularly in a problem with this degree of complexity.

Suppose P is a 2*N+1 x N+1 matrix of values in which the values in each k-th row are probability values. P(k,m) would be the probability of having m-1 unbroken pills remaining after k-1 days. Of course, the first row (zero days) would have a 1 in the last element and the rest would be zeros. You can create a for loop in which the outer loop starts with k = 2 and computes the next k-th row from the previous k-1 row. Within this loop, the k-th row can be calculated from the figures in the previous row. P(k,m) can be computed from P(k-1,m) and P(k-1,m+1). This is the heart of the computation. Think about it carefully. It depends on k as well as m.

After P is finally filled, it should be very simple to compute the expected value for each row in a single vector, and your task is done.