Calculate a probability exactly

Question

This task is about writing code to compute a probability exactly. The output should be a precise probability written as a fraction in its most reduced form. That is it should never output 4/8 but rather 1/2.

For some positive integer n, consider a uniformly random string of 1s and -1s of length n and call it A. Now concatenate to A its first value. That is A[1] = A[n+1] if indexing from 1. A now has length n+1. Now also consider a second random string of length n whose first n values are -1, 0, or 1 with probability 1/4,1/2, 1/4 each and call it B.

For example, consider n=3. Possible values for A and B could be A = [-1,1,1,-1] and B=[0,1,-1]. In this case the two inner products are 0 and 2.

Now consider the inner product of A[1,...,n] and B and the inner product of A[2,...,n+1] and B.

Your code must output the probability that both inner products are zero.

For n=1 this probability is clearly 1/2.

I don't mind how n is specified in the code but it should be very simple and obvious how to change it.

Languages and libraries

You can use any language and libraries you like. I would like to run your code so please include a full explanation for how to run/compile your code in linux if at all possible.

Answer 1

Pyth, 48 47 46 44 bytes

K,smlf!|Fms*Vd.>Tk2^,1_1Q^+0tM3Q^8Qj\//RiFKK

Try it online: Demonstration

The online version probably doesn't compute n=6. On my laptop (offline version) it takes about 45 seconds.

Brute force approach.

Explanation:

smlf!|Fms*Vd.>Tk2^,1_1Q^+0tM3Q   implicit: Q = input number
                          tM3    the list [-1, 0, 1]
                        +0       add zero, results in [0, -1, 0, 1]
                       ^     Q   all possible lists of length Q using these elements
 m                               map each list d (B in Lembik's notation) to:
                  ,1_1              the list [1, -1]
                 ^    Q             all possible lists of length Q
   f                                filter for lists T (A in Lembik's notation),
                                    which satisfy:
       m        2                      map each k in [0, 1] to:
        s*Vd.>Tk                          scalar-product d*(shifted T by k)
    !|F                                not or (True if both scalar-products are 0)      
  l                                 determine the length                
s                                add all possibilities at the end

K,...^8QQj\//RiFKK   
 ,...^8Q             the list [result of above, 8^Q]
K                    store it in K
              iFK    determine the gcd of the numbers in K
            /R   K   divide the numbers in K by the gcd
         j\/         join the two numbers by "/" and print

Answer 2

Mathematica, 159 100 87 86 85 bytes

n=3;1-Mean@Sign[##&@@Norm/@({1,0,0,-1}~t~n.Partition[#,2,1,1])&/@{1,-1}~(t=Tuples)~n]

To change n just change the variable definition at the beginning.

Since it's brute force it's fairly slow, but here are the first eight results:

n   P(n)
1   1/2
2   3/8
3   7/32
4   89/512
5   269/2048
6   903/8192
7   3035/32768
8   169801/2097152

The last one already took 231 seconds and the runtime is horribly exponential.

Explanation

As I said it's brute force. Essentially, I'm just enumerating all possible A and B, compute the two dot products for every possible pair and then find the fraction of pairs that yielded {0, 0}. Mathematica's combinatorics and linear algebra functions were quite helpful in golfing this:

{1,-1}~(t=Tuples)~n

This generates all n-tuples containing 1 or -1, i.e. all possible A. For n = 3 that is:

{{1, 1, 1}, 
 {1, 1, -1}, 
 {1, -1, 1}, 
 {1, -1, -1}, 
 {-1, 1, 1}, 
 {-1, 1, -1}, 
 {-1, -1, 1}, 
 {-1, -1, -1}}

To compute B we do almost the same:

{1,0,0,-1}~t~n

By repeating 0, we duplicate each tuple for each 0 it contains, thereby making 0 twice as likely as 1 or -1. Again using n = 3 as an example:

{{-1, -1, -1},
 {-1, -1, 0}, {-1, -1, 0},
 {-1, -1, 1},
 {-1, 0, -1}, {-1, 0, -1},
 {-1, 0, 0}, {-1, 0, 0}, {-1, 0, 0}, {-1, 0, 0},
 {-1, 0, 1}, {-1, 0, 1},
 {-1, 1, -1},
 {-1, 1, 0}, {-1, 1, 0},
 {-1, 1, 1},
 {0, -1, -1}, {0, -1, -1},
 {0, -1, 0}, {0, -1, 0}, {0, -1, 0}, {0, -1, 0},
 {0, -1, 1}, {0, -1, 1},
 {0, 0, -1}, {0, 0, -1}, {0, 0, -1}, {0, 0, -1},
 {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0},
 {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1},
 {0, 1, -1}, {0, 1, -1},
 {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
 {0, 1, 1}, {0, 1, 1},
 {1, -1, -1},
 {1, -1, 0}, {1, -1, 0},
 {1, -1, 1},
 {1, 0, -1}, {1, 0, -1},
 {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
 {1, 0, 1}, {1, 0, 1},
 {1, 1, -1},
 {1, 1, 0}, {1, 1, 0},
 {1, 1, 1}}

Now, for each possible A, we want the dot product of each of those possible B, both with A[1 .. n] and A[2 .. n+1]. E.g. if our current A is {1, 1, -1}, we want the dot product with both {1, 1, -1} and with {1, -1, 1}. Since all our B are already conveniently the rows of a matrix, we want the two sublists of A as columns of another matrix, so that we can compute a simple dot product between them. But transposing {{1, 1, -1}, {1, -1, 1}} simply gives {{1, 1}, {1, -1}, {-1, 1}} which is just a list of all 2-element cyclic sublists of A. That's what this does:

Partition[#,2,1,1]

So we compute that and take the dot product with our list of B. Since we now get a nested list (since each possible A yields a separate vector), we flatten those with ##&@@.

To find out if a pair {x, y} is {0, 0} we compute Sign[Norm[{x,y}]] where Norm gives √(x²+y²). This gives 0 or 1.

Finally, since we now just want to know the fractions of 1s in a list of 0s and 1s all we need is the arithmetic mean of the list. However, this yields the probability of both at least one dot product being non-zero, so we subtract it from 1 to get the desired result.

Answer 3

Pyth - 65 55 bytes

Fixed bug with fraction reduction at cost of one byte.

Uses brute force approach, and can be golfed hugely, but just wanted to get something out there. Very slow

*F-KP/Jmms*Vked,thdPhd*makhk^,1_1Q^[1ZZ_1)Q,ZZ2/lJ^2/K2

It uses Cartesian products to generate both A and B, doing the variable probabilities by making 0 appear twice in the source list and then counts the ones that inner product to zero. The inner product is make easy by the Vectorization syntactic sugar. Simplifing the fraction was scaring me initially, but it was pretty easy with the Prime Factorization function and the realization that we only have to reduce by powers of 2.

Try it online here.

Answer 4

CJam, 58 57 54 51 46 bytes

WX]m*Zm*_{~.+2,@fm<\f.*::+0-!},,__~)&:T/'/@,T/

To run it, insert the desired integer between WX] and m*.

Thanks to @jimmy23013 for the bit magic and for golfing off 5 bytes!

Try it online in the CJam interpreter.

Idea

Most parts of these answer are straightforward, but it uses two neat tricks:

• Instead of pairing all vectors of {-1, 1}ⁿ with all vectors of {-1, 0, 1}ⁿ with the desired probabilities, it considers counts the number of triplets of vectors in {-1, 1}ⁿ that satisfy a certain condition.

  If we add the last two vectors of a triplet, the result will be a vector of {-2, 0, 2}ⁿ.

  Since (-1) + 1 = 0 = 1 + (-1), 0s will occur twice as often as -2s and 2s.

  Dividing each component by 2 would yield a vector of {-1, 0, 1}ⁿ with the desired probabilities.

  Since we are only interested if the scalar product is 0 or not, we can skip the division by 2.

• After counting all triplets satisfying the condition of the question and the total number of triplets, we have to reduce the resulting fraction.

  Instead of calculating the GCD of both numbers, since the denominator will always be a power of 2, it suffices to divide both numbers by the highest power of 2 that divides the numerator.

  To obtain the highest power of 2 that divides x, we can take the bitwise AND of x and ~x + 1.

  ~x reverses all bits of x, so all trailing 0s become 1s. By adding 1 to ~x, those 1s will turn back into 0s and the last 1 in ~x + 1 will match the last 1 in x.

  All other bits are either both 0 of distinct, so the bitwise AND return the integer consisting of the last 1 of x and all 0s that follow it. This is the highest power of 2 that divides x.

Code

WX]    e# Push the array [-1 1].
       e# Insert N here.
m*     e# Cartesian product: Push the array of all vectors of {-1,1}^N.
Zm*    e# Cartesian product: Push the array of all triplets of these vectors.
_      e# Copy the array.
{      e# Filter; for each triplet of vectors U, V and W in {-1,1}^N:
  ~    e#   Dump U, V and W on the stack.
  .+   e#   Compute X := V + W, a vector of {-2,0,2}^N, where each component is
       e#   zero with probability 1/2.
  2,@  e#   Push [0 1]. Rotate U on top of it.
  fm<  e#   Push [U U'], where U' is U rotated one dimension to the left.
  \f.* e#   Push [U*X and U'*X], where * denotes the vectorized product.
  ::+  e#   Add the components of both products.
  0-   e#   Remove zeroes.
       e#   Push the logical NOT of the array.
},     e#   If the array was empty, keep the triplet.
,      e# Push X, the length of the filtered array.
__~)&  e# Push X & ~X + 1.
:T     e# Save the result in T and divide X by T.
'/     e# Push a slash.
@,T/   e# Dividet he length of the unfiltered array by T.