Wizard
Joined: Oct 14, 2009
• Posts: 25741
September 27th, 2011 at 5:18:42 PM permalink
Here is a question asked of me recently. Somebody I know hit six single-line royals in 5,000 hands once. He estimates he has played 25,000,000 hands of video poker in his lifetime. What are the odds of having such a streak at least once in 25M hands?

Of course, we could answer it by simulation, but that always leaves me so ... unsatisfied.

If we don't mind doing matrix multiplication 25 million times over, we could follow the kind of solution for a coin flipping problem I addressed in my June 4, 2010 Ask the Wizard column.

Any thoughts how to solve this? ME (MathExtremist), miplet, CM (CrystalMath), here is your chance to shine!
“Extraordinary claims require extraordinary evidence.” -- Carl Sagan
buzzpaff
Joined: Mar 8, 2011
• Posts: 5328
September 27th, 2011 at 5:38:44 PM permalink
.0003425
Wizard
Joined: Oct 14, 2009
• Posts: 25741
September 27th, 2011 at 6:05:50 PM permalink
Quote: buzzpaff

.0003425

“Extraordinary claims require extraordinary evidence.” -- Carl Sagan
miplet
Joined: Dec 1, 2009
• Posts: 2075
September 27th, 2011 at 6:15:50 PM permalink
Quote: Wizard

I'm pretty sure his work involved typing a random number. :) I have no clue how do it.
“Man Babes” #AxelFabulous
buzzpaff
Joined: Mar 8, 2011
• Posts: 5328
September 27th, 2011 at 6:22:49 PM permalink
NO ! Does speed count? I was first, you know. And so far no one has challenged my answer.

No fair miplet, how did you know ? LOL
weaselman
Joined: Jul 11, 2010
• Posts: 2349
September 27th, 2011 at 7:23:11 PM permalink
at least n royals in N hands is
F(n,N) = p*F(n-1,N-1) + (1-p)*F(n,N-1),
where p is 0.00000154

At least one such occurrence in 25M trials is P=1-(1-F(n,N))^25,000,000

I calculate F(6,5000) as 2.867*10^-16, and P = 8.327*10^-9

Seems too simple ... What am I missing?
"When two people always agree one of them is unnecessary"
Wizard
Joined: Oct 14, 2009
• Posts: 25741
September 27th, 2011 at 8:12:31 PM permalink
Quote: weaselman

where p is 0.00000154...Seems too simple ... What am I missing?

0.00000154 = 1 in 649,350. You seem to be looking at the royal probability on the deal. Let's use 1 in 40,000 for the probability on the draw.

That will get you closer, but there is still the issue that the 25,000,000 period of time overlap.
“Extraordinary claims require extraordinary evidence.” -- Carl Sagan
Wizard
Joined: Oct 14, 2009
• Posts: 25741
September 27th, 2011 at 8:32:47 PM permalink
Based on a quick estimate, my best guess at this time is 1 in 19,000.
“Extraordinary claims require extraordinary evidence.” -- Carl Sagan
Wizard
Joined: Oct 14, 2009
• Posts: 25741
September 28th, 2011 at 7:37:22 AM permalink
I ran a simulation overnight. Here are the results:

Wins=3
Losses=19125

Recall that each "loss" is 25M hands, so that is almost half a trillion hands.

So, based on this small sample size I get a probability of 1 in 6376. I'm still trying to find a mathematical solution, but no epiphany so far.
“Extraordinary claims require extraordinary evidence.” -- Carl Sagan
TIMSPEED
Joined: Aug 11, 2010
• Posts: 1246
September 28th, 2011 at 8:02:31 AM permalink
Heck, I hit 4 single-line royals in probably less than 1000 hands! (Different machines, but all inside of 4 days)
and I haven't hit one since : ( lol
Gambling calls to me...like this ~> http://www.youtube.com/watch?v=4Nap37mNSmQ
buzzpaff
Joined: Mar 8, 2011
• Posts: 5328
September 28th, 2011 at 9:06:35 AM permalink
R^{\alpha \beta} - {1 \over 2}R g^{\alpha \beta} + g^{\alpha \beta} \Lambda = \frac{8 \pi G}{c^4 \mu_0} ( F^{\alpha}{}^{\psi} F_{\psi}{}^{\beta} + {1 \over 4} g^{\alpha \beta} F_{\psi\tau} F^{\psi\tau}).

Or as stated before .0003425
CrystalMath
Joined: May 10, 2011
• Posts: 1909
September 28th, 2011 at 2:06:13 PM permalink
Quote: Wizard

Here is a question asked of me recently. Somebody I know hit six single-line royals in 5,000 hands once. He estimates he has played 25,000,000 hands of video poker in his lifetime. What are the odds of having such a streak at least once in 25M hands?

Of course, we could answer it by simulation, but that always leaves me so ... unsatisfied.

If we don't mind doing matrix multiplication 250 million times over, we could follow the kind of solution for a coin flipping problem I addressed in my June 4, 2010 Ask the Wizard column.

Any thoughts how to solve this? ME (MathExtremist), miplet, CM (CrystalMath), here is your chance to shine!

Wouldn't this be matrix multiplication 25 million times over and not 250 million?

If it is matrix multiplication, it actually won't be too bad, even if it is 25M or 250M. The tough part will be making the transition matrix.

Once we have that, we can calculate T^2. We can multiply that result by itself to get T^4. Multiply that result to get T^8, and so on until we get to T^(16,777,216). At this point, we've only done 24 multiplications. If we have stored all of the previous results, then we multiply the following matrices:

T^(16,777,216)
T^(4,194,304)
T^(2,097,152)
T^(1,048,576)
T^(524,288)
T^(262,144)
T^(65,536)
T^(16,384)
T^(8,192)
T^(4,096)
T^(2,048)
T^(64)

This will give us T^25,000,000

All together, this is 35 Matrix Multiplications.

Any thoughts on how to setup the transition matrix? Would we need to have 35,000 states (0 to 6 royals * 5000 previous hands)? The coin flipping was much easier because the streak had to occur back to back and not within the last 5,000 trials.
I heart Crystal Math.
CrystalMath
Joined: May 10, 2011
• Posts: 1909
September 28th, 2011 at 2:11:36 PM permalink
Or maybe we just need a 7 by 7 matrix? Still thinking...
I heart Crystal Math.
weaselman
Joined: Jul 11, 2010
• Posts: 2349
September 28th, 2011 at 2:26:40 PM permalink
The probability that the "winning" series start at hand k is P(k) = (1-P(k-1))*p*F(5,4999)
(F(n,N) is defined in my earlier post - it is the probability of hitting at least n royals in N hands).
Then P(1) = pF, P(2) = (1-pF)*pF ... P(k) = (1-pF)^(k-1)*pF
Summing these probabilities for 0 < k <= N, get: P = 1 - (1-pF)^N, where N = 24995000 - the last hand where the win could occur.
I get 0.000142794205019281 or 1 in 7003.08531333589 using this formula and 1/40000 for p.
"When two people always agree one of them is unnecessary"
CrystalMath
Joined: May 10, 2011
• Posts: 1909
September 28th, 2011 at 3:13:34 PM permalink
Here is my matrix solution. Please comment:

I started out with the initial state matrix populated with the binomial distribution of 0, 1, 2, 3, 4, 5, or 6+ royals in the first 5000 hands.

S0 = :
 0.8825 0.11031 0.0068935 0.00028712 8.9673e-06 2.2401e-07 4.7468e-09

For the transition matrix, this what I came up with:
 0.999975 2.5e-05 0 0 0 0 0 0.000199995 0.999775 2.4995e-05 0 0 0 0 0 0.00039999 0.999575 2.499e-05 0 0 0 0 0 0.000599985 0.999375 2.4985e-05 0 0 0 0 0 0.00079998 0.999175 2.498e-05 0 0 0 0 0 0.000999975 0.998975 2.4975e-05 0 0 0 0 0 0 1

Then I calculated T^(25,000,000 - 5,000):
 0.88238 0.1103 0.0068925 0.00028707 8.959e-06 2.1835e-07 0.00013625 0.88238 0.1103 0.0068925 0.00028707 8.959e-06 2.1835e-07 0.00013647 0.88237 0.1103 0.0068925 0.00028707 8.959e-06 2.1835e-07 0.00013843 0.88235 0.1103 0.0068923 0.00028706 8.9587e-06 2.1834e-07 0.00017008 0.88167 0.11021 0.006887 0.00028684 8.9519e-06 2.1817e-07 0.00093024 0.86019 0.10753 0.0067192 0.00027985 8.7338e-06 2.1286e-07 0.025275 0 0 0 0 0 0 1

Multiplying this by S0, I get:

 0.88238 0.1103 0.0068925 0.00028707 8.959e-06 2.1835e-07 0.00013632

With a final answer of: 0.000136318, or 1 in 7335.793909

I've got to run by son to soccer practice now, so details on the transition matrix will have to come later.
I heart Crystal Math.
Wizard
Joined: Oct 14, 2009
• Posts: 25741
September 28th, 2011 at 5:10:20 PM permalink
You, sir, are a GENIUS!

Outstanding work! That is the post of the month, easily.

I just spent the last hour confirming everything, and I agree exactly.

Thank you, very much, for your help.
“Extraordinary claims require extraordinary evidence.” -- Carl Sagan
dwheatley
Joined: Nov 16, 2009
• Posts: 1246
September 28th, 2011 at 6:39:17 PM permalink
I'm going to need the transition matrix explained to me. I don't immediately see how a one-step transition matrix can track how long into a streak we are. I was worried you would have to track how many hands ago each royal was, so you would know when 5000 hands have gone by and you have to drop one off the streak. This would require a lot of states (which I see was CM's initial instinct).
Wisdom is the quality that keeps you out of situations where you would otherwise need it
buzzpaff
Joined: Mar 8, 2011
• Posts: 5328
September 28th, 2011 at 7:47:54 PM permalink
No offense Crystal, but I stand by my answer. And I was first !
And my first four digits 0.000 agree with yours.
I hope you did not copy !
CrystalMath
Joined: May 10, 2011
• Posts: 1909
September 28th, 2011 at 7:59:33 PM permalink
Thanks Wizard, and you are very welcome.

So, here was the epiphany:

If p = the probability of a royal (assumed to be 1/40000) and you currently have n-royals in the last 5000 spins, where n is 1 to 5, then the following events can occur with one more game played:

drop a royal from the beginning hit a new royal new royal count probability
yes yes n n/5000 * p
yes no n-1 n/5000 * (1-p)
no yes n+1 (1-n/5000) * p
no no n (1-n/5000) * (1-p)

You see that the probability of dropping a royal is simply n/5000, since the royals will have an equal distribution in the last 5000.
I heart Crystal Math.
CrystalMath
Joined: May 10, 2011
• Posts: 1909
September 28th, 2011 at 8:03:11 PM permalink
Quote: buzzpaff

No offense Crystal, but I stand by my answer. And I was first !
And my first four digits 0.000 agree with yours.
I hope you did not copy !

You were amazingly close, well within an order of magnitude.

I heart Crystal Math.
CrystalMath
Joined: May 10, 2011
• Posts: 1909
September 28th, 2011 at 8:08:20 PM permalink
I didn't add the other cases:

If you have 0 royals, you can only gain a new one or not. The formulas hold true, though since the likelihood of dropping a royal is 0.
If you have 6 royals, you can change, but we don't care. Instead, we keep the player in that state so we can accumulate all hits in the 25 million games.
I heart Crystal Math.
buzzpaff
Joined: Mar 8, 2011
• Posts: 5328
September 28th, 2011 at 8:12:33 PM permalink
I did not see anything in the question Ruling out a WAG!! LOL
Wizard
Joined: Oct 14, 2009
• Posts: 25741
September 29th, 2011 at 9:35:48 AM permalink
Here is a preview of my explanation in my next column. Any comments?
“Extraordinary claims require extraordinary evidence.” -- Carl Sagan
miplet
Joined: Dec 1, 2009
• Posts: 2075
September 29th, 2011 at 10:07:21 AM permalink
Quote: A t W

Step 2 ......have 0, 1, 2, 3, 4, or 5 wilds

should say royals instead of wilds.
“Man Babes” #AxelFabulous
CrystalMath
Joined: May 10, 2011
• Posts: 1909
September 29th, 2011 at 10:07:44 AM permalink
Great write up. I like how you really spell it out and make it clear.

Just a few things:

1. In step 2, you wrote "0, 1, 2, 3, 4, or 5 wilds." Should be "royals."

2. In step 4, you list S-0, but you shortened the numbers, so some appear as zero. I would put the exact information shown in step 1, but I see you need to save space, so maybe scientific notation.
I heart Crystal Math.
Wizard
Joined: Oct 14, 2009
• Posts: 25741
September 29th, 2011 at 10:40:13 AM permalink
“Extraordinary claims require extraordinary evidence.” -- Carl Sagan
CrystalMath
Joined: May 10, 2011
• Posts: 1909
September 29th, 2011 at 10:50:18 AM permalink
I just read Franks question in the post: "What are the odd?"
I heart Crystal Math.
dwheatley
Joined: Nov 16, 2009
• Posts: 1246
September 29th, 2011 at 10:51:56 AM permalink
Your question (in bold) says "what are the odd".
Wisdom is the quality that keeps you out of situations where you would otherwise need it
DorothyGale
Joined: Nov 23, 2009
• Posts: 639
October 15th, 2011 at 3:41:42 PM permalink
I just now found this thread that relates to "Ask the Wizard! No. 277," where Mr. W. answers this question:
Quote:

I once hit six royals in single-line video poker within 5,000 hands. In my lifetime I have played about 25 million hands. What are the odds? Frank

The answer he arrives at is a "1 in 7,336 chance."

So, I decided to double check this by simulation ... how to simulate?

I looked up the frequency of a Royal in JOB 9/6, and based on that I used 1-in-40391 for a royal ... I then created a program to run 250M rounds of VP and see if 6 royals occur consecutively in 5000 hands at any point ... if so, you have a success, otherwise a failure ... I then put that program and a loop to run another 250M and then another and then another ... so, in the end I ran 218542 samples of 250M hands (well, slightly less than 250M, whenever I got 6 royals in 5k, I called that a success and exited) which amounted to a total of approximately 5.5T rounds of video poker.

Here are my results:

Number of trials of 250M hands: 218542
Number of trials where 6 royals occurred within 5000 hands: 32
Frequency of trials with 6 royals in 5000 hands: 1-in-6829.

According to Mr. W.'s work, the expected number of successes in 218542 is 29.8 with standard deviation 5.5, so my result is easily in the first standard deviation from expectation.

I obviously need to run a lot more trials of 250M to have this converge more closely ... but honestly, my computer is overheating ... I created 10 simultaneous background processes on my Ubuntu machine and let it run until I got tired of it.

Here is the code if you want to check it ... I didn't really check it very carefully, it could be wrong ... I'd be interested in a substantial increase in speed and efficiency if anyone can share that ...

#include <stdio.h>#include <stdlib.h>#include <time.h>main() {  srand(time(NULL));  int currentCount = 0;  int ctr, p, n;  int x[5000] = {};  int t = 0;  int success = 0;  for (n = 1; n <= 1000000; n++) {      currentCount = 0;    for (ctr = 0; ctr < 5000; ctr++)      x[ctr] = 0;    for (ctr = 0; ctr < 25000000; ctr++) {      p = rand()%40391;      if (p == 0 && x[ctr%5000] == 1)         currentCount = currentCount; // dummy      else      if (p == 0 && x[ctr%5000] == 0) {        currentCount++;        x[ctr%5000] = 1;      }      else      if (p > 0 && x[ctr%5000] == 1) {        currentCount--;        x[ctr%5000] = 0;      }      else      if (p > 0 && x[ctr%5000] == 0) {        currentCount = currentCount; // dummy      }      else        exit(1); // should never occur      if (currentCount >= 6) {        success++;        break;      }    }        printf("%6d, %6d, %1.6f\n", n, success, (double)success/n);    fflush(stdout);  }}
I believe my simulation confirms the result (not that it needed confirmation).

Good fun, and very nice solution Mr. W. and Crystal Math ...

--Ms. D.
"Who would have thought a good little girl like you could destroy my beautiful wickedness!"
weaselman
Joined: Jul 11, 2010
• Posts: 2349
October 15th, 2011 at 3:54:32 PM permalink
Quote: DorothyGale

I'd be interested in a substantial increase in speed and efficiency if anyone can share that ...

1. Rewrite the if/elseif statement to get rid of redundant operations. Particularly, avoid calculating the (ctr%5000) more than once, and get rid of "dummy" branches. These are relatively expensive things that you do over and over in a gazillion of iterations.
Something like
int idx = ctr%5000;if (p){   if(x[idx]) { x[idx] = 0; currentCount--; }}else{   if(!x[idx]) { x[idx]=1; currentCount++; }}

or, maybe, even:

int *xp = x + ctr%5000;if((*xp == 0) == (p == 0)) { *xp = !p; currentCount += p ? -1 : 1; }

should do the same stuff a bit faster.

2. Better yet, lose the x[5000] alltogether. Instead have a 6 element array, and fill it cyclically with the numbers of hands in which there is a hit. When filing the ith element, check the (i-5)th (wrapping around at 0, of course), if difference is less than 5000, count it as success, otherwise, keep going. That will save you some time on initialization in the beginning of each loop, and also %5 is faster than %5000.
If you insist on having x[5000], at least replace the loop filling it with zeroes with a call to memset().

FWIW, here is a sim I wrote in perl back when I was trying to confirm my solution to the original thread (I was using a recurrent function approach to avoid tedious matrix multiplication):
#!/usr/bin/perl$p = 1/40000;my ($total,$wins) = (0,0);while(1){$total++;   my $strikes = []; printf("%10d: %10f\n",$total, $wins/$total) unless $total % 100000; for (1 .. 25000000) { next unless rand() <=$p;      push(@$strikes,$_);      shift(@$strikes) if($_ - $strikes->[0] > 5000); ($wins++,last) if (@\$strikes >= 6);   }}
"When two people always agree one of them is unnecessary"
DorothyGale
Joined: Nov 23, 2009
• Posts: 639
October 15th, 2011 at 4:29:56 PM permalink
Mr. Weasel,

I'm not worried about the bits of rand(), I don't think it matters in this sim ... maybe it does, but hmmm .... I did #1 and the run time for 10 trials dropped from 7873416 microseconds to 6533979 microseconds ... I've got to say, though, that I was really looking for a new idea, like your perl ... can you time that for me for 10 trials? I ran it from the command line and it doesn't seem to work ...

Nothing came of these ... (yes I have perl installed)

My prompt> perl Yourprogram
My prompt> ./Yourprogram

--Ms. D.
"Who would have thought a good little girl like you could destroy my beautiful wickedness!"
CrystalMath
Joined: May 10, 2011
• Posts: 1909
October 15th, 2011 at 8:06:27 PM permalink
I just re-did the calculation with your video poker odds of 1:40,391 instead of 1:40,000 which we originally used. The probability under this scenario is 1:7765 .
I heart Crystal Math.
weaselman
Joined: Jul 11, 2010
• Posts: 2349
October 16th, 2011 at 9:06:16 AM permalink
Quote: DorothyGale

I did #1 and the run time for 10 trials dropped from 7873416 microseconds to 6533979 microseconds ... I've got to say, though, that I was really looking for a new idea, like your perl ... can you time that for me for 10 trials?
I ran it from the command line and it doesn't seem to work ...

That's because it only prints output every 100000 iterations :) You have to be really patient :)
It does not make very much sense to time the perl implementation, it is going to be at least an order of magnitude slower than "C". Interpreted languages are very slow doing massive explicit iterations.
Here is an equivalent program I wrote in "C". It takes about 0.38 seconds per iteration on my (fairly old) laptop.
#include<stdio.h>#include <stdlib.h>#include <time.h>struct sb { int x[6], start, end; };void shift(struct sb *b) { b->start = b->start == 5 ? 0 : b->start+1; }int first(struct sb *b) { return b->x[b->start]; }void push(struct sb *b, int v){    b->x[b->end]=v;    b->end = b->end == 5 ? 0 : b->end+1;}int sz(struct sb *b){    int s = b->end - b->start;    return s < 0 ? s+5 : s;}main(){    struct sb strikes;    int total = 0, wins = 0, hand = 0;    clock_t start = clock();    float elapsed = 0;    srandom(time(0));    while(1)    {        strikes.start = strikes.end = 0;        for(hand = 0; hand < 25000000; hand++)        {          if(random()%40000) { continue; }          push(&strikes, hand);          if(hand - first(&strikes) > 5000) { shift(&strikes); }          if(sz(&strikes) >= 6)          {              wins++;              break;          }        }        total++;        elapsed = ((float) clock() - start)/CLOCKS_PER_SEC;        printf("%6d, %10.6f, time: total %10.6f, avg: %10.6f\n", total, (float) wins/total, elapsed, elapsed/total);    }}
"When two people always agree one of them is unnecessary"
DorothyGale
Joined: Nov 23, 2009
• Posts: 639
October 16th, 2011 at 9:16:47 AM permalink
Mr. Weasel,

I timed your program for 10 samples of 25M -- your time in u-seconds is: 4691362. Pretty good 8-)

I compiled it with gcc and ran it under Ubuntu Linux 11.10 on my quad processor thing.

--Ms. D.
"Who would have thought a good little girl like you could destroy my beautiful wickedness!"
weaselman
Joined: Jul 11, 2010
• Posts: 2349
October 16th, 2011 at 10:20:45 AM permalink
How do you compile (what flags), and what exactly do you use to measure the time? What time does the program output itself report?

The reason I ask is that my 0.38 sec. number from a debuggable executable on a pretty old mac laptop, is already way lower than yours. I also ran an optimized (-O3) exec on a dual core 1.6G desktop, and got 0.22 seconds per iteration - more than twice faster than what you clocked. I wonder where the discrepancy comes from ...

    96,   0.000000,          0, time: total  21.230000, avg:   0.221145    97,   0.000000,          0, time: total  21.590000, avg:   0.222577    98,   0.000000,          0, time: total  21.830000, avg:   0.222755    99,   0.000000,          0, time: total  22.059999, avg:   0.225102   100,   0.000000,          0, time: total  22.290001, avg:   0.222900real    0m22.318suser    0m22.280ssys     0m0.020s
"When two people always agree one of them is unnecessary"
DorothyGale
Joined: Nov 23, 2009
• Posts: 639
October 16th, 2011 at 12:03:25 PM permalink
Quote: weaselman

How do you compile (what flags), and what exactly do you use to measure the time? What time does the program output itself report?

With -O3 it is faster, Total u-seconds: 4077101. Well that's on my dedicated Ubuntu machine that's 3 years old. On my Cygwin shell on my Windows 7 machine that's 2 years old I get 2043000. Why would Ubuntu be so slow?
"Who would have thought a good little girl like you could destroy my beautiful wickedness!"
weaselman
Joined: Jul 11, 2010
• Posts: 2349
October 16th, 2011 at 3:24:18 PM permalink
Hard to tell. Certainly, it is not about ubuntu vs. windoze, must be that either something else is running on your linux box, or a hardware difference. What cpu does each of the boxes have?
"When two people always agree one of them is unnecessary"
mustangsally
Joined: Mar 29, 2011
• Posts: 2463
September 18th, 2014 at 9:43:58 AM permalink
Quote: CrystalMath

Here is my matrix solution. Please comment:

sweet as candy

did you raise the matrix in Excel or use another program?

I set mine up in Excel too, but then made the matrix to use in R (save time)
I found some add-ins to raise a matrix in Excel but have not looked at them yet.

Quote: CrystalMath

With a final answer of: 0.000136318, or 1 in 7335.793909

nice
0.0001363178
1 in 7335.798

and R code (my homework for this week. please do not laugh at me, with me is okay)
I got lazy and used your matrix for rows 3:7 as I have more homework to do than calculate probabilities in Excel right now
the matrix and vector were taken right from Excel, so that is nice too
require (expm)p <- 1/40000   #I was playing with the matrix to see what it can dooooooooq <- 1-pA <- matrix(c(q,p,0,0,0,0,0,(1/5000)*q,(4999/5000)*q + (1/5000)*p,(4999/5000)*p,0,0,0,0,0,0.00039999,0.99957502,0.00002499,0,0,0,0,0,0.000599985,0.99937503,0.000024985,0,0,0,0,0,0.00079998,0.99917504,0.00002498,0,0,0,0,0,0.000999975,0.99897505,0.000024975,0,0,0,0,0,0,1), byrow=T, nrow=7,dimnames = list(c(0:5,"6+"),c(0:5,"6+ royals")))Asystem.time(B <-A%^%(25000000-5000))Bv0<-c(0.88249552366128,0.110314698325118,6.89346204814203E-03,2.87119872301923E-04,8.96733669524682E-06,2.24009670889038E-07,4.66231783332393E-09)C <- v0%*%BCp1 <- C[1,7]p11/p1

results
> require (expm)Loading required package: expmLoading required package: MatrixLoading required package: latticeAttaching package: ‘expm’The following object is masked from ‘package:Matrix’:    expm> > p <- 1/40000> q <- 1-p> > A <- matrix(c(q,p,0,0,0,0,0,+ (1/5000)*q,(4999/5000)*q + (1/5000)*p,(4999/5000)*p,0,0,0,0,+ 0,0.00039999,0.99957502,0.00002499,0,0,0,+ 0,0,0.000599985,0.99937503,0.000024985,0,0,+ 0,0,0,0.00079998,0.99917504,0.00002498,0,+ 0,0,0,0,0.000999975,0.99897505,0.000024975,+ 0,0,0,0,0,0,1), byrow=T, nrow=7,+ dimnames = list(c(0:5,"6+"),c(0:5,"6+ royals")))> A             0          1           2          3           4          50  0.999975000 0.00002500 0.000000000 0.00000000 0.000000000 0.000000001  0.000199995 0.99977501 0.000024995 0.00000000 0.000000000 0.000000002  0.000000000 0.00039999 0.999575020 0.00002499 0.000000000 0.000000003  0.000000000 0.00000000 0.000599985 0.99937503 0.000024985 0.000000004  0.000000000 0.00000000 0.000000000 0.00079998 0.999175040 0.000024985  0.000000000 0.00000000 0.000000000 0.00000000 0.000999975 0.998975056+ 0.000000000 0.00000000 0.000000000 0.00000000 0.000000000 0.00000000    6+ royals0  0.0000e+001  0.0000e+002  0.0000e+003  0.0000e+004  0.0000e+005  2.4975e-056+ 1.0000e+00> system.time(B <-A%^%(25000000-5000))   user  system elapsed       0       0       0 > B           0         1           2            3            4            50  0.8823753 0.1102997 0.006892508 0.0002870711 8.958995e-06 2.183479e-071  0.8823751 0.1102996 0.006892507 0.0002870710 8.958994e-06 2.183479e-072  0.8823734 0.1102994 0.006892493 0.0002870704 8.958976e-06 2.183475e-073  0.8823455 0.1102959 0.006892275 0.0002870613 8.958692e-06 2.183405e-074  0.8816746 0.1102121 0.006887035 0.0002868431 8.951881e-06 2.181745e-075  0.8601908 0.1075265 0.006719218 0.0002798536 8.733750e-06 2.128583e-076+ 0.0000000 0.0000000 0.000000000 0.0000000000 0.000000e+00 0.000000e+00      6+ royals0  0.00013625161  0.00013646972  0.00013843323  0.00017007954  0.00093024465  0.02527461446+ 1.0000000000> v0<-c(0.88249552366128,0.110314698325118,6.89346204814203E-03,2.87119872301923E-04,8.96733669524682E-06,2.24009670889038E-07,4.66231783332393E-09)> C <- v0%*%B> C             0         1           2           3            4            5[1,] 0.8823753 0.1102996 0.006892508 0.000287071 8.958995e-06 2.183479e-07        6+ royals[1,] 0.0001363178> C[1,7]   6+ royals 0.0001363178 > 1/C[1,7]6+ royals  7335.798
not all online R programs do not load the expm package
this one works
http://www.r-fiddle.org/#/
Loading required package: expmLoading required package: MatrixAttaching package: 'expm'The following object is masked from 'package:Matrix':    expm             0          1           2          3           4          50  0.999975000 0.00002500 0.000000000 0.00000000 0.000000000 0.000000001  0.000199995 0.99977501 0.000024995 0.00000000 0.000000000 0.000000002  0.000000000 0.00039999 0.999575020 0.00002499 0.000000000 0.000000003  0.000000000 0.00000000 0.000599985 0.99937503 0.000024985 0.000000004  0.000000000 0.00000000 0.000000000 0.00079998 0.999175040 0.000024985  0.000000000 0.00000000 0.000000000 0.00000000 0.000999975 0.998975056+ 0.000000000 0.00000000 0.000000000 0.00000000 0.000000000 0.00000000    6+ royals0  0.0000e+001  0.0000e+002  0.0000e+003  0.0000e+004  0.0000e+005  2.4975e-056+ 1.0000e+00           0         1           2            3            4            50  0.8823753 0.1102997 0.006892508 0.0002870711 8.958995e-06 2.183479e-071  0.8823751 0.1102996 0.006892507 0.0002870710 8.958994e-06 2.183479e-072  0.8823734 0.1102994 0.006892493 0.0002870704 8.958976e-06 2.183475e-073  0.8823455 0.1102959 0.006892275 0.0002870613 8.958692e-06 2.183405e-074  0.8816746 0.1102121 0.006887035 0.0002868431 8.951881e-06 2.181745e-075  0.8601908 0.1075265 0.006719218 0.0002798536 8.733750e-06 2.128583e-076+ 0.0000000 0.0000000 0.000000000 0.0000000000 0.000000e+00 0.000000e+00      6+ royals0  0.00013625161  0.00013646972  0.00013843323  0.00017007954  0.00093024465  0.02527461446+ 1.0000000000             0         1           2           3            4            5[1,] 0.8823753 0.1102996 0.006892508 0.000287071 8.958995e-06 2.183479e-07        6+ royals[1,] 0.0001363178   6+ royals 0.0001363178 6+ royals  7335.798 >

still Sally says stupid that Excel, R and many other programs can not do this properly and/or easily
A%^%(25000000-5000)
did you see the system time in R?

this was kind of fun

now another vacation!
Sally
I Heart Vi Hart