# Compute the Wilson score interval

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. nS and nF are the number of successes and failures, respectively, and n is the total number of trials (equivalent to nS + nF). 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)1.

Given non-negative integers nS and nF, 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


1. 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.

# 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!)

• Can Mathematica use arbitrary encodings? The square root symbol is 1 byte in a lot of single-byte encodings. – Mego Jan 3 '17 at 3:50
• 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. – Greg Martin Jan 3 '17 at 3:51
• Oh, it requires a command line switch (or some function call)? Gross. – Mego Jan 3 '17 at 3:52
• Mathematica is a wonderful juxtaposition of awesome and gross :D – Greg Martin Jan 3 '17 at 3:52

# 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.

• 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. – Mego Jan 3 '17 at 8:00

# 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


# 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.

# R, 58535149 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  • I like it when R is competitive with golfing languages... – J.Doe Oct 2 '18 at 19:38 # 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]. # 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! # 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²  • – Adám Oct 30 '18 at 22:24 • @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: ⊣×⊢÷+ -> ×÷+ – ngn Nov 1 '18 at 21:33 • @ngn Right, but actually, this is allowed by double meta consensus: (1) and (2). – Adám Nov 1 '18 at 22:28 # 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. # 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


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


# 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.

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


# 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)


# 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
½      sqrt[ ((s×f)÷(s+f)+W)×© ]   (call this R)

×Ø-            [   -R,      +R  ]

÷           Divide both by:


## 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 Ø-.

# 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