Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
Related: Simulating Theatre puzzle –  Mr.Wizard Oct 22 '13 at 17:52

1 Answer 1

up vote 4 down vote accepted

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 $$

share|improve this answer
    
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

Your Answer

 
discard

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.