6

I have discrete L-length line filled randomly with 2-length segments so its cover line with gaps 0 or 1. So we can describe cover configuration as sequence of gaps such as {0,0,1,0,1}. How I can calculate the probability of any cover? I tried to do it using weighted graphs but it works only for small values of L < 30. Increase the line length on 1 leads to 10-times increase in calculation time. Here is my slow code: 

l = 2;
detach = Function[{cov, l}, s = Last[cov]; 
Table[Append[cov, Delete[ReplacePart[s, i -> l + Total[s[[i ;; i + 1]]]], i + 1]], {i, 1, Length[s] - 1}]];
getProbability = Function[{cov, l},
   cc = NestWhile[Flatten[(detach[#, l] &) /@ #, 1] &, {cov}, Length@Last@Last@# != 1 &];
   Total[Map[(Apply[Times, #]) &, MapThread[(p = Total[Select[#1, # >= l &] - l + 1]; 
   If[ p === 0, 1, 1/p]) &, {cc,ConstantArray[ConstantArray[l, Length[First[First[cc]]]], Length[cc]]}, 2]]]
 ];
Table[getProbability[{ConstantArray[0, i]}, l] // AbsoluteTiming, {i, l, 10}]

Output: {{0.031540, 1}, {0.000216, 2/3}, {0.000547, 7/15}, {0.002399, 34/105}, {0.013129, 638/2835}, {0.089803, 4876/31185}, {0.694457, 220217/2027025}, {6.223379, 6885458/91216125}, {63.263433, 569311642/10854718875}}

Here is example graph:

enter image description here

Maybe it is possible to solve this problem using recurrence sequence or generating functions?

VividD
  • 3,660
  • 4
  • 26
  • 42
Филипп Цветков
  • 3,461
  • 1
  • 15
  • 34

1 Answers1

6

There is an elegant resursive formula which can be efficiently implemented with memoization

ClearAll[f];

f[gaps_List, l_: 2] := f[gaps, l] = Sum[f[gaps[[;; n]], l] f[gaps[[n + 1 ;;]], l], 
    {n, Length[gaps] - 1}]/(Total[l + gaps] - 2 l + 1);
f[{_}] = f[{_}, _] = 1;

Here Total[l + gaps] - 2 l + 1 is L - l + 1 (the number of avaliable positions on the empty line)

l = 2;
f[ConstantArray[0, #], l] & /@ Range[l, 10] // AbsoluteTiming

{0.000684, {1, 2/3, 7/15, 34/105, 638/2835, 4876/31185, 220217/2027025, 
   6885458/91216125, 569311642/10854718875}}

f[RandomInteger[l - 1, 30], l] // AbsoluteTiming

{0.062516, 1031/2604060900000}

It is really fast!

Proof

Let us consider final filling. For example ($L=7,\ l=2$)

enter image description here

One of these segments was chosen first. For example, we choose the segment, which is marked by gray.

enter image description here

Probability of these choice is $1/(L-l+1)=1/6$ (it doesn't depend on the position). Then we need to fill the empty space to obtain the final filling

enter image description here

It splits to two independent problems

enter image description here

In the gap notation this step can be written as

$$ (0,1,0,0) \to (0,1), (0,0) $$

However, there are another possibilities:

$$ (0,1,0,0) \to (0), (1,0,0)\\ (0,1,0,0) \to (0,1), (0,0)\\ (0,1,0,0) \to (0,1,0), (0) $$

As a result

$$ p(0,1,0,0) = \frac{1}{6}\bigl(p(0)p(1,0,0)+p(0,1)p(0,0)+p(0,0,1)p(0)\bigr) $$

In general:

$$ p(g_1,g_2,\ldots,g_n) = \frac{1}{L-l+1}\sum_{i=1}^{n-1}p(g_1,\ldots,g_i)p(g_{i+1},\ldots,g_n) $$

Here length $L$ depends on gaps as $$ L = \sum_{i=1}^{n}(l+g_i)-l $$

ybeltukov
  • 43,673
  • 5
  • 108
  • 212
  • For the simpletons (me) would you please tell how this works? – Mr.Wizard Oct 22 '13 at 17:52
  • Could you give me the reference on this elegant recursive formula in the book or journal publecation? – Филипп Цветков Oct 22 '13 at 18:21
  • @ybeltukov It is really cool! Is it works correctly for any distributions of the gaps? – Филипп Цветков Oct 22 '13 at 20:02
  • @Mr.Wizard Surely! I added the proof to my post. – ybeltukov Oct 22 '13 at 21:39
  • @ФилиппЦветков I don't know any references, I have derived this formula by myself for half an hour. I checked it with your function and results always was the same. – ybeltukov Oct 22 '13 at 21:44
  • I am wondering if I don't understand the question. As I read it this is like the "Theatre puzzle" I linked to (below the question), but looking for a specific filling rather than a proportion. By my method the specific filling you show (http://i.stack.imgur.com/8XuIC.png) should happen with a probability of ~0.187 yet I do not get a value similar to this from your function (though I may be misapplying it). – Mr.Wizard Oct 22 '13 at 22:00
  • @Mr.Wizard I get f[{0, 1, 0, 0}] = 3/16 = 0.1875. May be you use wrong gaps? – ybeltukov Oct 22 '13 at 22:04
  • You're right, I failed to understand the gap notation. – Mr.Wizard Oct 22 '13 at 22:07
  • Now that I can see that it is working: this is really nice. +1 By the way, if you have spare time you should try ProjectEuler.net if you have not done so already; you'd do very well at it! – Mr.Wizard Oct 22 '13 at 22:10
  • @ybeltukov It is really easy! Shame on me. I should think better before asking. – Филипп Цветков Oct 23 '13 at 10:04
  • @Mr.Wizard I alredy know about ProjectEuler.net (from your profile!). I used to think that I need to write programs in ordinary languages (I find it boring). Thanks for pointing me to understand it better! Spare time is a problem for be, but there are integersting problems. – ybeltukov Oct 24 '13 at 21:18