Compute the Wilson score interval

Question

The Wilson score interval is a confidence interval of the probability of success, based on the proportion of successes in a set of Bernoulli trials (a Bernoulli trial is a trial in which exactly two outcomes are possible: success or failure). The interval is given by the following formula:

The two values given by the formula are the upper and lower bounds of the interval. n_S and n_F are the number of successes and failures, respectively, and n is the total number of trials (equivalent to n_S + n_F). z is a parameter dependent on the level of confidence desired. For the purposes of this challenge, z = 1.96 will be used (corresponding to a 95% confidence interval)¹.

Given non-negative integers n_S and n_F, output the bounds of the Wilson score interval.

Rules

The outputs must be as accurate as possible to the true values, within the limits of your language's floating-point implementation, ignoring any potential issues due to floating-point arithmetic inaccuracies. If your language is capable of arbitrary-precision arithmetic, it must be at least as precise as IEEE 754 double-precision arithmetic.
The inputs will be within the representable range for your language's native integer type, and the outputs will be within the representable range for your language's native floating-point type.
n will always be positive.
The order of the outputs does not matter.

Test Cases

Format: n_s, n_f => lower, upper

0, 1 => 0.0, 0.7934567085261071
1, 0 => 0.20654329147389294, 1.0
1, 1 => 0.09452865480086611, 0.905471345199134
1, 10 => 0.016231752262825982, 0.3773646254862038
10, 1 => 0.6226353745137962, 0.9837682477371741
10, 90 => 0.05522854161313612, 0.1743673043676654
90, 10 => 0.8256326956323345, 0.9447714583868639
25, 75 => 0.17545094003724265, 0.3430464637007583
75, 25 => 0.6569535362992417, 0.8245490599627573
50, 50 => 0.40382982859014716, 0.5961701714098528
0, 100 => 0.0, 0.03699480747600191
100, 0 => 0.9630051925239981, 1.0

The z value is the 1-α/2th quantile of the standard normal distribution, where α is the significance level. If you want a 95% confidence interval, your significance level is α=0.05, and the z value is 1.96.

Related: Fastest Gun in the West Problem. I was going to make it a challenge, but I guess you beat me to it. :/ — mbomb007, Commented Feb 22, 2017 at 16:42

Greg Martin · Accepted Answer · 2017-01-03 03:49:44Z

6

Mathematica, 48 bytes (UTF-8 encoding)

({-1,1}√((s=1.4^4)##/+##+s^2/4)+#+s/2)/(s+##)&

Unnamed function taking two arguments in the order n_s, n_f and returning an ordered pair of real numbers. The three-byte symbol √, representing the square-root function, is U-221A.

Uses the fact that preceding ## by a number results in the product of the two arguments, while +## results in their sum. Also uses the fact that products and sums automatically thread over lists, so that {-1,1}√(...) implements the ± in the formula. Defining the constant s = z^2 instead of z itself also saved a couple of bytes. (Mostly I'm just proud of saving a byte by noticing that 1.4^4 is exactly 1.96^2!)

answered Jan 3, 2017 at 3:49

Greg Martin

16.2k4 gold badges21 silver badges72 bronze badges

\$\begingroup\$ Can Mathematica use arbitrary encodings? The square root symbol is 1 byte in a lot of single-byte encodings. \$\endgroup\$
– user45941
Commented Jan 3, 2017 at 3:50
\$\begingroup\$ It can indeed use many encodings—for example, Mac OS Roman, which has the property you mention. My understanding, though, is that I would need to include the bytes necessary to switch to a non-default encoding, which in this case is more than the 2 "wasted" bytes. \$\endgroup\$
– Greg Martin
Commented Jan 3, 2017 at 3:51
\$\begingroup\$ Oh, it requires a command line switch (or some function call)? Gross. \$\endgroup\$
– user45941
Commented Jan 3, 2017 at 3:52
4

\$\begingroup\$ Mathematica is a wonderful juxtaposition of awesome and gross :D \$\endgroup\$
– Greg Martin
Commented Jan 3, 2017 at 3:52

Add a comment |

Sean · Accepted Answer · 2017-01-03 07:27:17Z

3

Perl 6, 66 bytes

->\s,\f,\z=1.96 {(s+(-z|z)*sqrt(s*f/(s+f)+z*z/4)+z*z/2)/(s+f+z*z)}

This function actually returns an or-junction of the lower and upper bounds; for example, if called with the arguments 100 and 0, it returns:

any(0.963005192523998, 1)

It's a non-traditional output format to say the least, but no particular format was specified, and both of the required values are present.

answered Jan 3, 2017 at 7:27

Sean

8,03611 silver badges15 bronze badges

\$\begingroup\$ This seems fine to me. The lack of a specific output format was intentional - requiring a specific format for output gives advantages to some languages, and unnecessarily complicates the challenge. As long as both output values are present in some usable form, it's acceptable. \$\endgroup\$
– user45941
Commented Jan 3, 2017 at 8:00

Add a comment |

Emigna · Accepted Answer · 2017-01-03 10:41:22Z

05AB1E, 34 bytes

"1.96"Dn©4/¹P¹O/+t*D()®;+¹0è+®¹O+/

Input is of the form [n_s, n_f]
Output is of the form [upper, lower]

Try it online!

Explanation

"1.96"                             # push 1.96
      Dn©                          # duplicate, square, store a copy in register
         4/                        # divide by 4
           ¹P¹O/                   # product of input divided by sum of input
                +                  # add this to (z^2)/4
                 t*                # sqrt and multiply with z
                   D()             # wrap in list together with a negated copy
                      ®;+          # add (z^2)/2
                         ¹0è+      # add n_s
                             ®¹O+/ # divide by z^2+n

Draco18s no longer trusts SE · Accepted Answer · 2018-10-02 00:00:52Z

Runic Enchantments, 105 bytes

#StillBetterThanJava

/:2,:2,i:3s:3s:i:3s*3s+::6s4s,+'qA{*:Z4s3TRr4s{++}\
\p2:,C1Ä'<> yyyyyyyyyyyyyyyyyyy'Ä1C,2p:2,//@S ',+ /

Try it online!

Input is of the form n_s n_f
Output is of the form lower upper and has a trailing space

AH GOD this one's a mess. Here's the unwrapped version:

>'Ä1C,:2p:2,:2,i:3s:3s:i:3s*3s+::6s4s,+'qA{*:Z4s3TRr4s{++}+,' S@
                  > yyyyyyyyyyyyyyyyyyy'Ä1C,2p:2,//

All those ys are to slow down the second IP so that it arrives at the Transfer point at the right time (i.e. second). This shoves the top 3 items of the one pointer over to the other (the setup to this action is depicted below). 'Ä1C, generates z by dividing character 196 by 100 (dup, square, dup, div 2, dup, div 2...). The everything-else is just a bunch of math and stack manipulation to shove future values down the stack until they're needed. For the most part, they end up in the right order and its only until r4s{++} that we have to reverse the stack and rotate the whole thing to get the values we want next to each other next to each other.

There's probably room for improvement, but its complex enough that I can't see it. Heck, had inadvertently read "z" instead of "n" in the original formula at one point and fixing that was rough.

I had to pull out notecards and simulate the stacks in order to make sure it was correct:

Every single one has a value on both ends due to how many variables there were (eg. I'd have one with S and one with F, I'd pop them both, flip one around and add the S+F that was on the other end to the top of the stack). You can see one of the sqrt(...) cards has an S on the lower edge.

Robert S. · Accepted Answer · 2018-10-02 19:04:57Z

3

R, 58 53 51 49 41 bytes

-15 bytes thanks to J.Doe. -2 bytes thanks to Giuseppe.

function(x,y)prop.test(t(c(x,y)),cor=F)$c

edited Oct 2, 2018 at 19:04

answered Oct 1, 2018 at 21:02

Robert S.

1,32311 silver badges17 bronze badges

1

\$\begingroup\$ I like it when R is competitive with golfing languages... \$\endgroup\$
– J.Doe
Commented Oct 2, 2018 at 19:38

Add a comment |

Luis Mendo · Accepted Answer · 2017-01-03 02:26:46Z

2

MATL, 38 37 bytes

pGs/1.96XJU4/XI+X^J*t_hIE+G1)+GsJU+/S

Input is an array of two numbers, in any of these formats: [25 75], [25, 75], [25; 75].

Try it online! or verify all test cases.

edited Jan 3, 2017 at 2:26

answered Jan 3, 2017 at 2:20

Luis Mendo

105k9 gold badges135 silver badges372 bronze badges

Add a comment |

flawr · Accepted Answer · 2017-01-03 11:14:44Z

2

Haskell, 70 69 68 67 bytes

s#f|y<-1.96^2=[(s+y/2+k*sqrt(s*f*y/(s+f)+y^2/4))/(s+f+y)|k<-[-1,1]]

Try it online!

edited Jan 3, 2017 at 11:14

answered Jan 3, 2017 at 10:55

flawr

43.9k7 gold badges104 silver badges249 bronze badges

Add a comment |

Community · Accepted Answer · 2020-06-17 09:04:33Z

2

APL (Dyalog Unicode), 50 bytes

{(+/⍺⍵z)÷⍨(⍺+z÷2)(-,+).5*⍨z×(⍺×⍵÷⍺+⍵)+4÷⍨z←3.8416}

Try it online!

Infix Dfn, taking \$⍺←n_s\$ and \$⍵←n_f\$.

Thanks to H.PWiz and dzaima for helping out.

How:

                                        z←3.8416  ⍝ Set z=1.96²
                                     4÷⍨          ⍝ Divide it by 4
                                    +             ⍝ plus
                           (⍺×⍵÷⍺+⍵)             ⍝ (nf×ns)÷n
                         z×                       ⍝ ×z²
                     .5*⍨                         ⍝ to the power .5 (square root)
                (-,+)                             ⍝ ±
         (⍺+z÷2)                                  ⍝ ns+(z²/2)
(+/⍺⍵z)÷⍨                                        ⍝ all divided by nf+ns+z²

edited Jun 17, 2020 at 9:04

CommunityBot

1

answered Oct 5, 2018 at 13:51

J. Sallé

3,4539 silver badges24 bronze badges

\$\begingroup\$ 40: z←1.9208⋄(z+⊣(-,+).5*⍨z×z+2×⊣×⊢÷+)÷z+z++ \$\endgroup\$
– Adám
Commented Oct 30, 2018 at 22:24
\$\begingroup\$ @Adám That's neither an expression, nor a complete program, but you could make it an expression by initialising z at its rightmost usage: ...÷z+(z←1.908)++ for the same byte count. Also: ⊣×⊢÷+ -> ×÷+ \$\endgroup\$
– ngn
Commented Nov 1, 2018 at 21:33
1

\$\begingroup\$ @ngn Right, but actually, this is allowed by double meta consensus: (1) and (2). \$\endgroup\$
– Adám
Commented Nov 1, 2018 at 22:28

Add a comment |

orlp · Accepted Answer · 2017-01-03 13:04:18Z

1

Python, 79 67 bytes

lambda s,f,z=3.8416:2j**.5*(s-(-z*(f*s/(f+s)+z/4))**.5+z/2)/(f+s+z)

Output is a complex integer with the interval stored as the real/imaginary part.

edited Jan 3, 2017 at 13:04

answered Jan 3, 2017 at 9:33

orlp

39.2k7 gold badges81 silver badges195 bronze badges

Add a comment |

Community · Accepted Answer · 2020-06-17 09:04:33Z

1

dc, 71 bytes

16k?dsa2*?sb1.96 2^dso+dlalb4**lalb+/lo+vlov*dsd+lalblo++2*dsx/rld-lx/f

Takes both inputs on two separate lines upon invocation, and outputs on two separate lines with the upper bound on the bottom and the lower bound on top.

For example:

bash-4.4$ dc -e '16k?dsa2*?sb1.96 2^dso+dlalb4**lalb+/lo+vlov*dsd+lalblo++2*dsx/rld-lx/f'
10                # Input n_s
90                # Input n_f
.0552285416131361 # Output lower bound
.1743673043676654 # Output upper bound

edited Jun 17, 2020 at 9:04

CommunityBot

1

answered Jan 6, 2017 at 5:18

R. Kap

4,9082 gold badges15 silver badges37 bronze badges

Add a comment |

rnso · Accepted Answer · 2017-01-09 03:58:16Z

1

Racket 134 bytes

(let*((n(+ s f))(z 1.96)(h(* z z))(p(/ 1(+ n h)))(q(+ s(/ h 2)))(r(* z(sqrt(+(/(* s f) n)(/ h 4))))))(values(* p(- q r))(* p(+ q r))))

Ungolfed:

(define (g s f)
  (let* ((n (+ s f))
         (z 1.96)
         (zq (* z z))
         (p (/ 1 (+ n zq)))
         (q (+ s (/ zq 2)))
         (r (* z (sqrt (+ (/(* s f) n) (/ zq 4))))))
    (values (* p (- q r)) (* p (+ q r)))))

Testing:

(g 1 10)

Output:

0.016231752262825982
0.3773646254862038

answered Jan 9, 2017 at 3:58

rnso

1,70710 silver badges8 bronze badges

Add a comment |

Community · Accepted Answer · 2020-06-17 09:04:33Z

1

Java 7, 130 bytes

Golfed:

double[]w(int s,int f){double n=s+f,z=1.96,x=z*z,p=s+x/2,d=z*Math.sqrt(s*f/n+x/4),m=1/(n+x);return new double[]{m*(p-d),m*(p+d)};}

Ungolfed:

double[] w(int s, int f)
{
    double n = s + f, z = 1.96, x = z * z, p = s + x / 2, d = z * Math.sqrt(s * f / n + x / 4), m = 1 / (n + x);
    return new double[]
    { m * (p - d), m * (p + d) };
}

Try it online

Returns an array of type double of length 2, can probably be golfed more.

edited Jun 17, 2020 at 9:04

CommunityBot

1

answered Jan 9, 2017 at 13:47

peech

3091 silver badge4 bronze badges

Add a comment |

Sok · Accepted Answer · 2018-10-02 12:12:01Z

><> with `-v` flag, 100 bytes

:{:@:}*@+:@,24a,a,-:@:*:&4,+:\$~*{&:&2,+}:{:@@-{&+:&,nao+&,n;
,}:{::*@@-:0$0(?$-1a,:*:*:*(?\}:{:@,+2

Expects the input to be present on the stack at execution start, in the order n_s, n_f. Try it online!

What a stupid language to attempt this in...

As ><> lacks an exponent or root operator, the square root is calculated in the second line of code using the Babylonian method, to an accuracy of 1e-8 - for every example I've tried, this is accurate to at least 10 decimal places. If this isn't precise enough, the bounds can be tightened by adding more :* in the second line, shuffling things around to keep the mirrors in line.

Output is in the following form:

<lower bound>
<upper bound>

Community · Accepted Answer · 2020-06-17 09:04:33Z

Pyth, 38 bytes

J*K1.96Kmc++dhQcJ2+sQJ_B*K@+cJ4c*FQsQ2

Input is as a list of values, [n_s, n_f]. Output is [upper, lower] Try it online here, or verify all the test cases at once here.

J*K1.96Kmc++dhQcJ2+sQJ_B*K@+cJ4c*FQsQ2   Implicit: Q=eval(input())

  K1.96                                  Set variable K=1.96 (z)
J*K    K                                 Set variable J=K*K (z^2)
                                *FQ      Product of input pair
                               c   sQ    Divide the above by the sum of the input pair
                            cJ4          J / 4
                           +             Add the two previous results
                          @          2   Take the square root of the above
                        *K               Multiply by K
                      _B                 Pair the above with itself, negated
        m                                Map each in the above, as d, using:
             hQ                            First value of input (i.e. n_s)
               cJ2                         J / 2
          ++d                              Sum the above two with d
         c                                 Divided by...
                  +sQJ                     ... (J + sum of input)

lynn · Accepted Answer · 2018-10-02 15:20:10Z

Jelly, 30 bytes

×÷++1.96²©HH¤×®½×Ø-+®H¤+³÷++®ɗ

Try it online!

Explanation

                   Inputs: s and f
×÷+                (s×f)÷(s+f)
    1.96²©HH¤      (© ← 1.96²)÷4      (call this W)

   +                         (s×f)÷(s+f)+W
             ×®             ((s×f)÷(s+f)+W)×©
               ½      sqrt[ ((s×f)÷(s+f)+W)×© ]   (call this R)

                ×Ø-            [   -R,      +R  ]
                   +®H¤        [©/2-R,   ©/2+R  ]
                       +³      [©/2-R+s, ©/2+R+s]

                         ÷           Divide both by:
                          ++®ɗ       (s+f)+©

Note

Some of these features are newer than the challenge. I believe around the time this challenge was posted, ++®¶×÷++1.96²©HH¤×®½×-,1+®H¤+³÷ç was valid Jelly (32 bytes), lacking ɗ and Ø-.

Community · Accepted Answer · 2020-06-17 09:04:33Z

1

APL(NARS), 49 chars, 98 bytes

{k←1.9208⋄(2+n÷k)÷⍨(1+⍺÷k)+¯1 1×√1+2×⍺×⍵÷k×n←⍺+⍵}

test

  f←{k←1.9208⋄(2+n÷k)÷⍨(1+⍺÷k)+¯1 1×√1+2×⍺×⍵÷k×n←⍺+⍵}
  25 f 75
0.17545094 0.3430464637 
  0 f 1
0 0.7934567085 
  1 f 0
0.2065432915 1

edited Jun 17, 2020 at 9:04

CommunityBot

1

answered Oct 4, 2018 at 10:01

user58988

Add a comment |

Stack Exchange Network

Compute the Wilson score interval

Rules

Test Cases

16 Answers 16

Mathematica, 48 bytes (UTF-8 encoding)

Perl 6, 66 bytes

05AB1E, 34 bytes

Runic Enchantments, 105 bytes

R, 58 53 51 49 41 bytes

MATL, 38 37 bytes

Haskell, 70 69 68 67 bytes

APL (Dyalog Unicode), 50 bytes

How:

Python, 79 67 bytes

dc, 71 bytes

Racket 134 bytes

Java 7, 130 bytes

><> with `-v` flag, 100 bytes

Pyth, 38 bytes

Jelly, 30 bytes

Explanation

Note

APL(NARS), 49 chars, 98 bytes

Your Answer

Hot Network Questions

Compute the Wilson score interval

Rules

Test Cases

16 Answers 16

Mathematica, 48 bytes (UTF-8 encoding)

Perl 6, 66 bytes

05AB1E, 34 bytes

Runic Enchantments, 105 bytes

R, 58 53 51 49 41 bytes

MATL, 38 37 bytes

Haskell, 70 69 68 67 bytes

APL (Dyalog Unicode), 50 bytes

How:

Python, 79 67 bytes

dc, 71 bytes

Racket 134 bytes

Java 7, 130 bytes

><> with -v flag, 100 bytes

Pyth, 38 bytes

Jelly, 30 bytes

Explanation

Note

APL(NARS), 49 chars, 98 bytes

Your Answer

Sign up or log in

Post as a guest

Related

Hot Network Questions

><> with `-v` flag, 100 bytes