<jupyter_start><jupyter_text>Sequential testing
A common issue with classical A/B-tests, especially when you want to be able to detect small differences, is that the sample size needed can be prohibitively large. In many cases it can take several weeks, months or even years to collect enough data to conclude a test. 

*   The lower number of errors we require, the larger sample size we need.
* The smaller the difference we want to detect, the larger sample size is required.

Sequential sampling works in a very non-traditional way; instead of a fixed sample size, you choose one item (or a few) at a time, and then test your hypothesis. You can either:

* Reject the null hypothesis (H0) in favor of the alternate hypothesis (H1) and stop,
* Keep the null hypothesis and stop,
* unable to reach either conclusion with current observation and continue sampling.

 Advantage of Sequential testing over classic a/b testing


*   optimize necessary observation (sample size)
*   reduce the likelihood of error
* gives a chance to finish experiments earlier without increasing the possibility of false results

**N.B: Unlike classical fixed sample-size tests, where significance is only checked after all samples have been collected, this test will continously check for significance at every new sample and stop the test as soon as a significant result is detected, while still guaranteeing the same type-1 and type-2 errors as the fixed-samplesize test.**

 Common sequential testing algorithms
- [Evan Miller Simple Sequential A/B Testing](https://www.evanmiller.org/sequential-ab-testing.html)
- [Sequential probability ratio testing (SPRT)](https://en.wikipedia.org/wiki/Sequential_probability_ratio_test)

The Evan Miller sequential procedure for one-sided test works as follows:
* choose a sample size \(N\). [Here](https://www.evanmiller.org/ab-testing/sequential.html) is the link to define N
* Assign subjects randomly to the treatment and control, with 50% probability each.
* Track the number of incoming successes from the treatment group. Call this number \(T\).
* Track the number of incoming successes from the control group. Call this number \(C\).
* If \(T-C\) reaches \(2\sqrt{N}\), stop the test. Declare the treatment to be the winner.
* If \(T+C\) reaches \(N\), stop the test. Declare no winner.

The two-sided test is essentially the same, but with an alternate ending:
* If \(T-C\) reaches \(2.25\sqrt{N}\), stop the test. Declare the treatment to be the winner.
* If \(C-T\) reaches \(2.25\sqrt{N}\), stop the test. Declare the control to be the winner.
* If \(T+C\) reaches \(N\), stop the test. Declare no winner.

**N.B: This test completely ignores the number of failures in each group, which makes it significantly easier to implement in low-conversion settings. However**
- **If we hit the threshold without having reached statistical proof, we cannot continue the experiment.**
- **Even if preliminary evidence says that one version is terrible, we will keep losing conversions until we hit an arbitrarly pre-agreed threshold.**

 Sequential probability ratio testing (SPRT)
SPRT is based on the likelihood ratio statistic.
 Variants of SPRT
- [Wald's(1947) sequential probability ratio test](http://en.wikipedia.org/wiki/Sequential_probability_ratio_test)
 * Wald's test uses one pair of observations at each stage and ignores tied observations and thus reduces a test of a single proportion.
- [Conditional sequential probability ratio test](http://www.jstor.org/stable/2346379?seq=1page_scan_tab_contents)
 * used to calculate tied observations b/n two proportions. So, determining statistical significance with two distributional streams of data can be used for conditional sprt.
- [Sequential Generalized Likelihood Ratio Test](https://www.auduno.com/2014/12/25/rapid-a-b-testing-with-sequential-analysis/)
  * used to measure tied observation. Has stoping rule estimation problem.

We focus on conditional SPRT for this challenge.

General steps of conditional SPRT

1. Calculate critical upper and lower decision boundaries
2. Perform cummlative sum of the observation
3. calculate test statistics(likelihood ration) for each of the observations 
4. calculate upper and lower limits for exposed group
5. apply stopping

 Stopping Rule

1. If the log probability ratio greater than or equal to the upper critical limit then the model reject the null hypothesis with the favor of alternative hypothesis (i.e. accept H1 (conclude that version two is better than version one)) and terminate the test.
2. If the log probability ratio less than or equal to the lower critical limit then accept the null hypothesis (i.e. conclude that there is no difference between the two groups) and terminate the test.
3. If neither critical limit is reached, conduct another trial and continue the test.

One-sided SPRT

![image.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAPsAAAClCAYAAACa9oshAAAUBElEQVR4Ae2ZC5Ikqa5Ee2m9tN75jKkzT7WXCgggRPxSmE25vo5wiMl57/76lSsVSAVSgVQgFUgFUoFU4HMV+O/Xr1/+nx416OmpLdVYP2svFzyJqUAq0FDAf2jeb7ROp47YY3q4bEwFnqqA//Bqfi2OLuRB4obESqh5enyd1miO+sRUIBXoUICPh1L11ba8+jW7p057W/Waa/VYXa5UIBXYUICPSJEWjalteXy1iYGag7MUa9WTA+HxPvHEVCAVqCigH43aVu59pdCc2lpT4/D16qut/bW43y/9VCAVqCjQ+oh6c6N1vfU2MrUgx/A+8cRUIBWoKFD6aDSGDUJT833c6okZsoiZr7b6Gldba+BLTAVSgUUK+I9v0TZJmwqkAmcqkB/6mern3qlAKpAKpAKpQCqQCqQCqUAqkAqkAqnAS4FV//f5FO+fP3+sL<KEY>)

Two-sided SPRT

![image.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAP4AAAChCAYAAADnTgJzAAATw0lEQVR4Ae2Y25LcOA5E5/9/ejfgrmOnYUCiSEqiVOBDJ64JMKmKCc9//9UpBUqBUqAUKAVKgVKgFCgFSoFSoBQoBUqBUqAUKAVKgVIgUOB/QaxCpUAp8GIF7EdfP/wXP3BdrRSIFKgffqRKxUqBlytQP/yXP3BdrxSIFKgffqRKxUqBlytQP/yXP3BdrxSIFKgffqRKxUqBlytQP/yXP3BdrxSIFKgffqRKxUqBlyuw0g+fXcBR6Y/yUA+Ozrd+uAzrlALLKMCHefdCfg/vH92PfkOO2sRA6jOf+BFUTrWPcFRtKXCKAit8kHs7kAdNCGxFBNKY2XrI7cWiPL0Ran1k0xPlKlYKXK7ACh/k3g5R3sfUN9uOxj6h36A5tX8XiOHze760/jZ9z+9EGV+lwDLfwQqL+B3wDe3gf9ww1lKT9ftefEM7+B9316cO9P3EC9+nAG+9h7ffnAXvXCTaQWNqs6ePed/qZsU8z57PjtkOmi/7eQrw/iPYdWsb6E8U8zWRz/JR7sqY30N9tdnJx7xvdT5mfnaiWuq3csbn88zI4uQL11SAd5uN3bfdW6SHGM6e3tk97ALC732LE1OkHiSHv4fUg9RnPnHD6GgeO6qr2PUK8B69GG2ccUW1h2OzyeE7vMjNDXfufefsm2V/zHjeqBdbL5rxt/YfqmMYTeb3Hs/Vy3N135173zn7ap1XnccbjODI3bK5I5y7vQyl0PvEW3Ckt4W/akqBXgX4Nnuxd27Wl+2R1U+P2wL+RDFfE/lcJspVrBQ4UwG+vV48czflzvbTmkvsmYvAdcniNeSrFODbGsG7BMt2vmuf33P9Yr8TBw14DrZVeSnwSwG+n15cTcbsHqvtObwPFx0mKoJXKsD30YtPECW725K7s6wtZ3bvUZ5ejup7rgK8/wg+9fbZnZe9jy5sS+L3LDzS2zOveq5XgDfuxes3Pm9ipsF5Eycy2/J2QG//ZNv+IkRbdVWtqABv2Isr3mn2Tpk2s+ecymeXsAN6+yfb9hdB2qqr6i4FeKdevGvvO+dmWt2509DsmReCa2ihah5WgHcYweElXkCwpd9jr2eX4vgLEj+K8Bztq/pjCqDzKB6b+j3Vma6vUMAuZwf8uEOAYEMk1fxLAbQcxZKzTYFM57buB1VlFyXec5WR3p55T+5Bqxn4ZB3u3j3T/+69Tp2fXdriPQe+nt639aDFDHybNnffJ3uTu/e6ZL5d3g74cYcAQYdIHtTMfWfgg6792FWzd3rshXoWNxHsgB93CBB2iGShZu4zCxe62teskr3d1wjgL5oJQtzXt/gjvS38s2vYdwbO3q34xhTI3nSM9SXdmTgW7znw9fSe1cNOo3jWfsU7T4HsjedNeBGTimXXMr/3wNXb39PHzFHsmV09ayiQvf0a2y26hYpmK+L3rDvSuzUP3hHc4q/c8xTIvoXn3eSmjU1AO6C3f7Jtf3mMtuo/VfSN4B+2st6sQPaNvPnOp9zNhLQDevsn2/aXR/HVxEfQc5b/PQpk3833KHDCTWeKmnG1xE+42u9/tjB/1gz4DGedGVy61ww+vRvcGjvbZqbHs+d+Df8sYT2P+leLyWzmep/4UVQetY/yaP0Wj+VajufwfgtHVgNX6y4ZT0tcZ6nd0ls1NymgDxXZV63F7GweedDqsBWzfuLU4rfw0JOhcmFTi++xNa912Iqe13zNm33W8XPwz5pXvBMV8I+F73HiyJCKeWFSPmbN+x7vay12VONj6qttHPggvBlmdVkcnijvY97XXrOzPHU9CKfHHq7qaVTAxJ59eMCIl5zHqHY0xgx48A3t4H/cMBbVbNWT833qq231+CAcGWZ1Po5vaAf/44axqGarXnNHbWZ5PMpT9R0KmOh2wI87BDzkFgk1Hrd6jubg1j6NqU2Nj3mfOsMjOa1VW3l8XGdpnY/jR/0aUzvriWqoNdzLa21k0+8xqq3YiQr4B/B+z2g4Wnup99jav1UHJzXqqx3lLRbVbMVbeDwnPgiHouVajudQX224fMz71IF7eeo80ufR15V/oQL+MdTvWYP+o730eTzK4+szPuJaT0xR89iaxyZn6GOZn8WV66gNJ0i/93VPcoZbh7qtGnLUeiRfuIACew9+ZEUe+kiP1tIfodadYTPzDO4VOc+4L5weV7x/7ST/leLBekUZ7de5cHnUmpk2c2Zyrsw1677weFz57rVb8KPnAXvEGenN5sEZYdZT8fMViN7DYnUeokD0WFGs5Tp8DC21PTXwe+zhqp7jCnjd8Y8zVcftCtjj+RPFfE3kX/UhMCfCaK+KjSkQ6WyxOg9WYOajwnWlHMyM8Mo93jYr0tNidV6kgH/k3qvB09s/0sfsCEd4v6k30o7YN+lQdz2owCofCXtEePBKry+PNCL2+svXBecocMYH4zm9v7U5tRFu9X1DLtKE2Dfcv+44UYEzPhzP6f2W9enJsIXjDTXZ/S1e58sUiB49irXIwofVUttaAyc74dOPDxLPkLoIs56nx6O7WqzOFyqQfQzEeyQZ6c3mwWloB9/bkf+rIfkDT4ZJ22PC2b0sXufLFcg+jl5Z4Ovtj/qUE9vQDv7H/ccnvofwRLjXu1o+uoPF6pQCfykw86Pgo/trwKCjnNjsjM8I7xM/gnBEeITnytpoV2JX7lGzHqYAH4lizxXo7+nNejyn+mpbv/czzpY4XBG29F9RE+1msTqlwK4CMz8euHaHHijwnJnv4wdG7JbCHeFu8+SCaAdik0cV3ZsVsI9m1nn7B8j9IpylYcYTzbTYGw935W7eJ34ElUNt48A3jPxP+F3AZWfcCgFncK3OwV0jnLV7xE1s1owVefwdvd+zc8bh497vmfWIHi7qsWd5OHp6n9rDnSPsvVPEZbFvOXp/uzM+98cHiXskr6h8Go/sqHYkZjP06EziGvP11ExBPwi/h3ykt2feaj3cP8K9XaMeYnu9b8tzb0M7+N6O/F8NrsfXtfJpnXJkcT97q05z2KDnwZ+OOtDs3qM8vRxv6EOHCPV+UZ6Y1n2brRpgG9rB/7j/+Flc+zLb81OXoZ/VWufn7PGQn4q6bLZQ60C4WuqP1Hq+lt6WGs97hs8erXjGDk/jRCvbG9tQ/Y/7O48P0hf5mlPbatVXGx5fE/nU+n711c7qiZ+CtoAd0Ns/2ba/0WWyziO1GceT4tw3wifd44pd0YhZ6qttee9HPb5Oe9TeqtPcVk9rnXJgg8Zhx/uf8Bwwcjugt3+ybX/3FiWvyDwfY2IUJ7bVqzVZHTWKzJ2Fyt1iz5r7ZB504g6Z7+PUg+QVLYfv7S1/q0f7tuo053vMt0MN+AnPBwZ47JkER9Trc/ggPerv2Zq3fvUzW+u0RuPs0ovwRqicUZ6Y1pVdCpyiAB8b2Dtkq9/n8DO0Hcj5fYh71DpyEQ85kD7vE29F+iPc44h6iO31Vr4U6FKAD0yxh4j+qNfn8MGWHmp8Dz5odZmtOa3ROHNaEI4IW/qjmojLYnVKgWkKzPzI4IqW8zl8kB7192zNW7/6ma11WqNxdtlCeiPc6juSi7iJHeGp2lLgHwXsQ5p19j5K8oo2W32z9WiOODHzscGoJqvb64ELpD5Cas7CaKbF6pQCXQrM/Hj4OLsWWbiJe0V49drRDsSu3qXmPVgBPhqPPVeCo6d3tR7uEuEqu0a7WaxOKbCrwMyPB67doQsXcIcIV1072pXYqjvXXjcrYB/IrPPUj429I5ylzVU80R0sVqcU+EuBmR8FH91fAxZ12DXDRdduXiu7l8XrlAL//F9xPpgeaUZ6e+b19LBjhD18T+iJ7mqxOl+swMyPAq7V5GSvDFfb96x9svtbvM4XKsAHYVcf+QiUZwUZ2SfCFfa7c4dIE2J37lWzL1KAxza0g/9xD8FI76FBO8XsEeFO69elI42IfZ0Y33Rhe2Q7oLd/sm1/7/xgmB1h2/ZVFWlHrNR5mQL2sHZAb/9k2/7e8ZEwM8K2rasqUiDS02J1XqLAzAeG62xpmBPh2bO/kT/S2WJ1Hq6Af9je68DT27/XB7/Hvb7Kz1HA644/h71YLlXAHm/WOeNDgDPCWXsXz3EFovewWJ2HKDDzsfgYZlwdLo9HuEd69+bAvVfXmje+0cNO4Cgf/fBFO2pObXoLF1VAH0vtnnXp7+m1HvojPMoJB33eJ96DcBnOOPBFXK0zPIf3I+6WmPKoHfWS9xjVVuxmBfwj4fes1dtLn8eeHawHnqyfPKg9xAyjo/moxud9TZTX2NbMKGcx+vfyWoetmPUTpxY/Q+o8ZvUVv0EB/zj4Pasc7aXeY89s7YFPY2pHeR/zPv0Wt7OVp8bX+R588Ic5/5vVZXGYoryPeZ9exZYarTebHo++rvyLFfAPgt+zRksvNR575mU9cJPHN7SD/3HDWFSzVU/O96mvttXjg3BkmNX5OL6hHfyPG8aimq16zbXY8Hts6a2aixSwx+k5PGrUS85jVDsaY4byaExtanzM+9SBWd7H1VfbePBBuD325LVHbbh9zPvUGW7ltK7Vhs9ja3/VnaSAPUjP4SHpxfdI/kxkJjPUVzvKWyyqoXYr7/vUV1s5fNzPUT+zPYf6atPvY97P6ojPQGZ6nMFdHDsKeNHxd9rCNL0Zhk0nBv0ejCKOb0hMUfPepm4v7uvwQfq9T/wIwgHS632LE1OkXlHz2JqfZcPtcRZ/8TgFvNDm956Ii1gv51V9T9lzlh6r3pe9PM66d/GIAiayHcT+uIeB/h48PGxyAztPpl2W7gn3ZUePy4r6pMUQ1XbGNuw59Pte4iPoOcv/HgWy7+Z7FDjhpiaqHcTF/oke+6scRzrpG8Ej86r2uQpk38hzb3TT5iakHQTF/oke+6scxzq3q+Edwe0JlX2aAtm38LR73LavCmhL4PcsNNLbM0/3ZXYv9s6vvvsVyN78/s0W3wDhbE2ze4/y9HLM7mOnUZy9V/HNVyB74/mTXsIYCdZzNXh6eu/oYd8ZeMf+NTNXIHvTvOPLMjMFgustEnKfWfgWXZ50j+ztnnSHU3Y1YeyA3v7Jtv1F5Lbq51dx3xn4fDXWv0H2TutvfsKGJoYdL8onfAjgONT00mK0mIEvlei2a2VvcttCdw02IewgyMc9DKP9hwc+uAGtZuCDZbh99Uz/2xd70gKI+KSdV90VLUdx1futtlem82p7Tt2HSxup2b1HeXo5qm9fAXQexf1J31mR6foqNfSSdjH8nkuO9PbMq55YAd5hBGPm74pu6fd4JexydkBv/2Tb/iJUW3VV3aUA79SLd+1959xMqzt3GpptF7IDevsn2/YXcdqqq2pFBXjDXlzxTrN3yrSZPedUvpmXgOvUhYv8VgV44168dfnJwzMNJo85j85foHcSPL391fdsBXj/EXyqAtmdl7tPtijxnoVHenvmVc+zFOD76MUn3Da72zK7ZwsS71l0pLdnXvW8SwG+n15cTY3sHrfu6ZeasQycM7iMAz5whBcOw9ajPUf69vjh3aur/B8F0GwE/7Bda2U7X7uFmxYt5UqaXbiaGzYKPZf3N1r/SWmvt/8p/gS0zkLez/r24vAY1pmngOraY8/bZJsp22276+SsLtUziv6eXu3Z4yEPWi+2onJik8985aLGIxxb6HuUl76opmLnKIDmvTh7q2yP2XNSvmiBtHgnAddO2W56jyfK+5j3GdoSz2oyDl/vfe0zO8tTV3itArzHCI5snM0d4Ux7/bC08EACzgMtYannwTe0g/9xw1hLTdbve/EN7eB/3F2fOtD3Ey9cUwHeqxdbb5Xxt/Y31WVDiDeRuKKRXqWKeDSmNn0+tudnfRb3vT7m83s+s0BfT7zwmQrwnr0Y3TrjimoPxTJi4ofIPsUjvX6e51Jfbfp8TH21qTe0eHZ8j/pqW/+e72f4ep8v/z0K8NazcSmFuNyspeAD4fW+xYkp+vooR02EWm82h7j3iWstNYrUaazs71SAb2EEb1eO5e9Y5Ftn36F1zbxGAb7pPbxmm40pLLhRclrqW2efJmgRL6/And/8X+Iss8hfW5VTCpQCpypQP/xT5S3yUmBNBeqHv+a71FalwKkK1A//VHmLvBRYU4H64a/5LrVVKXCqAvXDP1XeIi8F1lSgfvhrvkttVQqcqkD98E+Vt8hLgTUVqB/+mu9SW5UCpypQP/xT5S3yUmBNBeqHv+a71FalwOkK2I+/TilQCpQCpUApUAqUAqVAKVAKlAIPVYD/nwGudA3951bLfi01V9xP91D7itk1oxTYVcB/lN7fJTixoGeXnp4zrrDKHmfcrThfoMDeB0oe1CsTi5A6cuZjg76GeITaH/X5GD7oOYkbas77mlM7q4OXWq1Tm7whR2PY5EDiiuQKS4FmBfQDMlsPOWLqq215fND3ZHHtPdrj6yOuKKa7mM0hDvq497fqNJfZxkcO9DPwwdY66gtLgU0F+KBAK8b2qDlIoxqt83l8rWnl2qpT3q06cob0eKSGuPezuHJu2Zrb4mKu1hPzfcQLS4FDCuiHpLYn8Tl8kHp8kLiiz+GD1KrfYtMH0gNaPLI15mvU36rTXGa3clkdR7ks5n3qCkuBTQX8h6O+2kaivtqaOxrXXhaFA/Rx36N1akd9Lb2eI/OzeMsMv9sWF7WeN/K1tuxSYFMBPjpQi4mBWzmrsUOtYhT/qf5TH/lwaL+v26rZqs1yytcyl3pDDjHzMzvLUa98Ea/vp6awFLhMAT7WywYuMGjmnf2PfCb3AlLVCm9V4Bs/1Nl3hg9867dS9yoFSoFS4LgC/wc6ZT3AKdtsRQAAAABJRU5ErkJggg==)

Next we will write python implementation for conditional SPRT<jupyter_code>#Importing Libraries
import sys
import os
import pandas as pd
from datetime import timedelta
import math
import numpy as np
import Pylance
#importing local scripts
#Adding scripts path
sys.path.append(os.path.abspath(os.path.join('..')))
#Importing data manipulator scripts
from scripts.data_manipulator import get_group, get_answered<jupyter_output><empty_output><jupyter_text>Importing Working Data<jupyter_code>#Reading Data
data= pd.read_csv('../data/AdSmartABdata.csv')
data.head(5)
data['date'] = pd.to_datetime(data['date'])
# # data['hour'] = pd.to_datetime(data['hour'], format='%H')
# data['hour'] = [f':{t}' for t in data['hour']]
# # data['hour'] = pd.to_datetime(data['hour'])
data.head()

#Function that builds a bernouli dataframe
def beronlli_constructor(engagement:pd.Series, yes_count:pd.Series):
    engagement = engagement.tolist()
    yes_count = yes_count.tolist()
    ret_list = []
    for i in range(len(engagement)):
        ret_list.extend(yes_count[i] * [1])
        ret_list.extend((engagement[i] - yes_count[i]) * [0])
        
    return ret_list
def transform_data(df):
    '''
        segment data into exposed and control groups
        consider that SmartAd runs the experment hourly, group data into hours. 
            Hint: create new column to hold date+hour and use df.column.map(lambda x:  pd.Timestamp(x,tz=None).strftime('%Y-%m-%d:%H'))
        create two dataframes with bernouli series 1 for posetive(yes) and 0 for negative(no)
        Hint: Given engagement(sum of yes and no until current observation as an array) and success (yes countas an array), the method generates random binomial distribution
            #Example
                engagement = np.array([5, 3, 3])
                yes = np.array([2, 0, 3])       
                Output is "[1] 1 0 1 0 0 0 0 0 1 1 1", showing a binary array of 5+3+3 values
                of which 2 of the first 5 are ones, 0 of the next 3 are ones, and all 3 of
                the last 3 are ones where position the ones is randomly distributed within each group.
    '''
    # Get answered quesionnaires only
    clean_data = get_answered(df=df, column1='yes', column2='no')

    #Segmenting data into exposed and control groups
    control_group = get_group(df=clean_data,from_column='experiment',classification='control')
    exposed_group = get_group(df=clean_data,from_column='experiment',classification='exposed')

    #Grouping Data using hours
    #For the control group
    control_group['date-time'] = control_group['date'] + control_group['hour'].apply(lambda x: timedelta(hours=x))
    #For the exposed group
    exposed_group['date-time'] = exposed_group['date'] + exposed_group['hour'].apply(lambda x: timedelta(hours=x))

    #Generating engagment and yes array based on date-time column grouping
    #control Group
    for_control_group = pd.DataFrame()
    for_control_group['engagement'] = control_group.groupby('date-time').agg({'yes':'count'}).iloc[:,0]
    for_control_group['yes_count'] = control_group.groupby('date-time').agg({'yes': 'sum'}).iloc[:,0]
    #Build bernouli series
    control_group_bernouli = beronlli_constructor(for_control_group['engagement'],for_control_group['yes_count'])

    #exposed Group
    for_exposed_group = pd.DataFrame()
    for_exposed_group['engagement'] = exposed_group.groupby('date-time').agg({'yes':'count'}).iloc[:,0]
    for_exposed_group['yes_count'] = exposed_group.groupby('date-time').agg({'yes': 'sum'}).iloc[:,0]
    #Build bernouli series
    exposed_group_bernouli = beronlli_constructor(for_exposed_group['engagement'],for_exposed_group['yes_count'])

    data_stream = (control_group_bernouli, exposed_group_bernouli)

    return data_stream
data_stream = transform_data(data)
#REFERENCE
# A Conditional Sequential Test for the Equality of Two Binomial Proportions
# <NAME>, Jr
# Journal of the Royal Statistical Society. Series C (Applied Statistics)
# Vol. 30, No. 2 (1981), pp. 109-115
class ConditionalSPRT:
    def __init__(self,exposed,control,odd_ratio,alpha=0.05,beta=0.10,stop=None):
        self.x = exposed
        self.y = control
        self.t1 = odd_ratio
        self.alpha = alpha
        self.beta = beta
        self.stop = stop

    def ConditionalSPRT(self,x,y,t1,alpha=0.05,beta=0.10,stop=None):
        """
        #
        # Meeker's SPRT for matched `x` (treatment) and `y` (control), 
        # both indicator responses, likelihood ratio t1, error rates alpha and beta,
        # and (optionally) truncation after trial stop.
        #
        # The return variable contains these elements:
        #(outcome,n, k,l,u,truncated,truncate_decision,x1,r,stats,limits)
        # * outcome:   "continue," "reject null," or "accept null".
        # * n: number observation used for the decsion
        # * k:     Index at which the outcome decision was made (or NA)
        # * l:     lower critical point
        # * u:     upper critical point
        # * truncate_decision: The approximate decision made after truncate point
        # * truncated: If the test was truncated, the value of `n.0`; NA otherwise
        # * x1:       Original data `x`, cumulative
        # * r:         Cumulative sum of x+y
        # * stats:     Series of cumulative sums of log probability ratios
        # * limits:    Two rows giving lower and upper critical limits, respectively
        #
        """
        if t1<=1:
            printLog('warning',"Odd ratio should exceed 1.")
        if (alpha >0.5) | (beta >0.5):
            printLog('warning',"Unrealistic values of alpha or beta were passed."
                    +" You should have good reason to use large alpha & beta values")
        if stop!=None:
            stop=math.floor(n0)

        def comb(n, k):
            return math.factorial(n) // math.factorial(k) // math.factorial(n - k)
        
        def lchoose(b, j):
            a=[]
            if (type(j) is list) | (isinstance(j,np.ndarray)==True):
                if len(j)<2:
                    j=j[0]
            if (type(j) is list) | (isinstance(j,np.ndarray)==True):
                for k in j:
                    n=b
                    if (0 <= k) & (k<= n):
                        a.append(math.log(comb(n,k)))
                    else:
                        a.append(0)
            else:
                n=b
                k=j
                if (0 <= k) & (k<= n):
                    a.append(math.log(comb(n,k)))
                else:
                    a.append(0)

            return np.array(a)

        def g(x,r,n,t1,t0=1):
            """
            #
            # Meeker's (1981) function `g`, the log probability ratio.
            # 
            """
            return -math.log(h(x,r,n,t1))+math.log(h(x,r,n,t0))

        def h(x,r,n,t=1):
            """
            #
            # Reciprocal of Meeker's (1981) function `h`: the conditional probability of 
            # `x` given `r` and `n`, when the odds ratio is `t`.
            #
            # `x` is his "x1", the number of positives in `n` control trials.
            # `r` is the total number of positives.
            # `n` is the number of (control, treatment) pairs.
            # `t` is the odds ratio.
            #
            """
            return f(r,n,t,offset=ftermlog(x,r,n,t))

        def f(r,n,t,offset=0):
            """#
            # Meeker's (1981) function exp(F(r,n,t)), proportional to the probability of 
            #  `r` (=x1+x2) in `n` paired trials with an odds ratio of `t`.
            #
            # This function does *not* vectorize over its arguments.
            #"""
            upper=max(0,r-n)
            lower=min(n,r)
            rng=list(range(upper,lower+1))
            return np.sum(fterm(rng,r,n,t,offset))

        def fterm(j,r,n,t,offset=0):
            ftlog=ftermlog(j,r,n,t,offset)
            return np.array([math.exp(ex) for ex in ftlog])

        def ftermlog(j,r,n,t,offset=0):
            """
            #
            # Up to an additive constant, the log probability that (x1, x1+x2) = (j, r) 
            # in `n` paired trials with odds ratio of `t`.
            #
            # `offset` is used to adjust the result to avoid under/overflow.
            #
            """
            xx=r-j
            lch=lchoose(n,j)
            lchdiff=lchoose(n,xx)
            lg=np.array(j)*math.log(t)
            lgsum=lch+lchdiff
            lgsum2=lgsum+lg
            lgdiff=lgsum2-offset

            return lgdiff

        def logf(r,n,t,offset=0):
            """
            #
            # A protected vesion of log(f), Meeker's function `F`.
            #
            """
            z=f(r,n,t,offset)
            if z>0:
                return math.log(z)
            else:
                return np.nan

        def clowerUpper(r,n,t1c,t0=1,alpha=0.05,beta=0.10):
            """
            #
            # Meeker's (1981) functions c_L(r,n) and c_U(r,n), the  critical values for x1.
            # 0 <= r <= 2n; t1 >= t0 > 0.
            #
            """
            offset=ftermlog(math.ceil(r/2),r,n,t1c)
            z=logf(r,n,t1c,logf(r,n,t0,offset)+offset)
            a=-math.log(alpha/(1-beta))
            b=math.log(beta/(1-alpha))
            lower=b
            upper=1+a
            return (np.array([lower,upper])+z)/math.log(t1c/t0)

        l=math.log(beta/(1-alpha))
        u=-math.log(alpha/(1-beta))
        sample_size=min(len(x),len(y))
        n=np.array(range(1,sample_size+1))

        if stop!=None:
            n=np.array([z for z in n if z<=stop])
        x1=np.cumsum(x[n-1])
        r=x1+np.cumsum(y[n-1])
        stats=np.array(list(map(g,x1, r, n, [t1]*len(x1)))) #recurcively calls g
        #
        # Perform the test by finding the first index, if any, at which `stats`
        # falls outside the open interval (l, u).
        #
        clu=list(map(clowerUpper,r,n,[t1]*len(r),[1]*len(r),[alpha]*len(r), [beta]*len(r)))
        limits=[]
        for v in clu:
            inArray=[]
            for vin in v:
                inArray.append(math.floor(vin))
            limits.append(np.array(inArray))
        limits=np.array(limits)

        k=np.where((stats>=u) | (stats<=l))
        cvalues=stats[k]
        if cvalues.shape[0]<1:
            k= np.nan
            outcome='Unable to conclude.Needs more sample.'
        else:
            k=np.min(k)
            if stats[k]>=u:
                outcome=f'Exposed group produced a statistically significant increase.'
            else:
                outcome='Their is no statistically significant difference between two test groups'
        if (stop!=None) & (k==np.nan):
        #
        # Truncate at trial stop, using Meeker's H0-conservative formula (2.2).
        # Leave k=NA to indicate the decision was made due to truncation.
        #
            c1=clowerUpper(r,stop,t1,alpha,beta)
            c1=math.floor(np.mean(c1)-0.5)
            if x1[n0]<=c1:
                truncate_decision='h0'
                outcome='Maximum Limit Decision. The aproximate decision point shows their is no statistically significant difference between two test groups'
            else:
                truncate_decision='h1'
                outcome=f'Maximum Limit Decision. The aproximate decision point shows exposed group produced a statistically significant increase.'
            truncated=stop
        else:
            truncate_decision='Non'
            truncated=np.nan
        return (outcome,n, k,l,u,truncated,truncate_decision,x1,r,stats,limits)
#Calculate odd Value
control_df = pd.DataFrame(data_stream[0])
exposed_df = pd.DataFrame(data_stream[1])
#counting values
p1_no_count, p1_yes_count = exposed_df.value_counts().values
p2_no_count, p2_yes_count = control_df.value_counts().values

p1 =  p1_yes_count / (p1_no_count + p1_yes_count)
p2 =  p2_yes_count / (p2_no_count + p2_yes_count)

odd_value = (p1/ (1 - p1)) / (p2/ (1 - p2))
odd_value
# Initialize the ConditionalSPRT class
cond_sprt = ConditionalSPRT(exposed=np.array(data_stream[1]),control=np.array(data_stream[0]),odd_ratio=odd_value)
result = cond_sprt.ConditionalSPRT(x=cond_sprt.x,y=cond_sprt.y,t1=cond_sprt.t1)<jupyter_output><empty_output><jupyter_text>The return variable contains these elements:
        (outcome,n, k,l,u,truncated,truncate_decision,x1,r,stats,limits)
         * outcome:   "continue," "reject null," or "accept null".
         * n: number observation used for the decsion
         * k:     Index at which the outcome decision was made (or NA)
         * l:     lower critical point
         * u:     upper critical point
         * truncate_decision: The approximate decision made after truncate point
         * truncated: If the test was truncated, the value of `n.0`; NA otherwise
         * x1:       Original data `x`, cumulative
         * r:         Cumulative sum of x+y
         * stats:     Series of cumulative sums of log probability ratios
         * limits:    Two rows giving lower and upper critical limits, respectively<jupyter_code>#Our final Decision
result[0]
#Plotting the result
limits = result[10].tolist()
lower_limit = [value[0] for value in limits]
upper_limit = [value[1] for value in limits]

x1 = result[7].tolist()
r = result[8]
import matplotlib.pyplot as plt
plt.plot(r, upper_limit,color = 'blue', linewidth=2, linestyle='dashed',label='Upper Bound')
plt.plot(r, lower_limit,color = 'red', linewidth=2, linestyle='dashed',label='Lower Bound')
plt.plot(r, x1,color = 'orange', linewidth=2, linestyle='solid',label='Cumulative Sum')
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08920528284601878, 'fraction_numerical': 0.061416902053591876, 'mean_word_length': 3.802749770852429, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 22, 'https://': 4, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '324299', 'license': 'MIT', 'n_tokens_mistral': 11931, 'n_tokens_neox': 11008, 'n_words': 1777, 'path': 'notebooks/sequential_testing.ipynb', 'repo_name': 'DePacifier/SmartAd-Campaign', 'size': '32619'}