Calculate the probability of being a squid coin with Bayes' theorem [python]

problem

I have two coins One coin has a 1/2 chance of appearing The other is a squid coin with 70% facing up Win if you can guess whether it is a squid coin

When I tried 50 more times now, it was 35 times.

By the way, can we call it a squid coin at this time?

For classical statistics

If it is not a squid coin, it should be in the top 25 times out of 50 times. When the difference in ratio is tested, Chi-square value is 4.16 p-value is 0.04 If you test with a risk rate of 5%, you can't say that it's a non-squid coin because it's about 4% even if it happens that "throwing a non-squid coin will make a table 35 times". Will be interpreted as

Not an ambiguous result such as "I can't say it" I want a probability of 1,0 such as whether it is a squid coin

In Bayesian statistics

P(X|Y) = \frac{P(Y|X)P(X)}{P(Y)}

in this case Binomial distribution

\begin{align}
 {}_n C_{x} p^x (1-p)^{n-x}
\end{align}

Likelihood P (Y | X) Treat as

Prior probability P (X) puts the probability of being a squid coin or a fair coin

This time, both Ikasama coins and fair coins will be selected equally. P(X)=0.5 Think of it as.

Possible events are Ikasama coin is selected and the table appears 35 times Fair coins are selected and the table appears 35 times

If this is expressed by a calculation formula

** Fair **

\begin{align}
 {}_{50} C_{35} (0.5)^{35} (1-0.5)^{50-35}=0.0019
\end{align}

** Ikasama **

\begin{align}
 {}_{50} C_{35} (0.7)^{35} (1-0.7)^{50-35}=0.1223
\end{align}

Will be. Multiply each value by the probability that Ikasama coin will be selected as the prior probability Since this alone cannot be considered as a probability, add both values to obtain a marginal distribution.

** Fair **

\frac{0.1223×0.5}{0.1223×0.5+0.0019×0.5}=0.016

** Ikasama **

\frac{0.0019×0.5}{0.1223×0.5+0.0019×0.5}=0.984
import math
import numpy as np
import random

def combinations_count(n, r):
    return math.factorial(n) // (math.factorial(n - r) * math.factorial(r))

loaded_prior=0.5
loaded_heads=0.7
loaded_tails=0.3

fair_prior=1-loaded_prior
fair_heads=0.5
fair_tails=0.5

n=50
x=35
n_x=n-x

def coin(n,x,n_x,lp,lh,lt,fp,fh,ft):
    loaded_likelihood=combinations_count(n, x)*(loaded_heads**x)*(loaded_tails**n_x)
    fair_likelihood=combinations_count(n, x)*(fair_heads**x)*(fair_tails**n_x)
    marginal = loaded_likelihood*loaded_prior + fair_likelihood*fair_prior
    load_prob=(loaded_likelihood*loaded_prior)/marginal
    fair_prob=(fair_likelihood*fair_prior)/marginal
    return load_prob, fair_prob

load_prob,fair_prob= coin(n=n,x=x,n_x=n_x,lp=loaded_prior,lh=loaded_heads,lt=loaded_tails,fp=fair_prior,fh=fair_heads,ft=fair_tails)
print("loaded coin "+str(round(load_prob,3))+"%, fair coin "+str(round(fair_prob,3))+"%")
loaded coin 0.984%, fair coin 0.016%

Calculated in this way, the probability of being a squid coin was 0.5 in advance, but increased to 0.98. By the way, if the table appears 33 times, it will exceed 90%.

try_time = 50

for o in range(try_time):
    n_x=try_time-o
    lo_p,fa_p=coin(n=try_time,x=o,n_x=n_x,
                          lp=loaded_prior,lh=loaded_heads,lt=loaded_tails,
                          fp=fair_prior,fh=fair_heads,ft=fair_tails)
    print("head "+str(o)+"  loaded coin "+str(round(lo_p,3))+"%, fair coin "+str(round(fa_p,3))+"%")
head 0  loaded coin 0.0%, fair coin 1.0%
head 1  loaded coin 0.0%, fair coin 1.0%
head 2  loaded coin 0.0%, fair coin 1.0%
head 3  loaded coin 0.0%, fair coin 1.0%
head 4  loaded coin 0.0%, fair coin 1.0%
head 5  loaded coin 0.0%, fair coin 1.0%
head 6  loaded coin 0.0%, fair coin 1.0%
head 7  loaded coin 0.0%, fair coin 1.0%
head 8  loaded coin 0.0%, fair coin 1.0%
head 9  loaded coin 0.0%, fair coin 1.0%
head 10  loaded coin 0.0%, fair coin 1.0%
head 11  loaded coin 0.0%, fair coin 1.0%
head 12  loaded coin 0.0%, fair coin 1.0%
head 13  loaded coin 0.0%, fair coin 1.0%
head 14  loaded coin 0.0%, fair coin 1.0%
head 15  loaded coin 0.0%, fair coin 1.0%
head 16  loaded coin 0.0%, fair coin 1.0%
head 17  loaded coin 0.0%, fair coin 1.0%
head 18  loaded coin 0.0%, fair coin 1.0%
head 19  loaded coin 0.0%, fair coin 1.0%
head 20  loaded coin 0.0%, fair coin 1.0%
head 21  loaded coin 0.0%, fair coin 1.0%
head 22  loaded coin 0.001%, fair coin 0.999%
head 23  loaded coin 0.002%, fair coin 0.998%
head 24  loaded coin 0.005%, fair coin 0.995%
head 25  loaded coin 0.013%, fair coin 0.987%
head 26  loaded coin 0.029%, fair coin 0.971%
head 27  loaded coin 0.065%, fair coin 0.935%
head 28  loaded coin 0.14%, fair coin 0.86%
head 29  loaded coin 0.275%, fair coin 0.725%
head 30  loaded coin 0.469%, fair coin 0.531%
head 31  loaded coin 0.674%, fair coin 0.326%
head 32  loaded coin 0.828%, fair coin 0.172%
head 33  loaded coin 0.918%, fair coin 0.082%
head 34  loaded coin 0.963%, fair coin 0.037%
head 35  loaded coin 0.984%, fair coin 0.016%
head 36  loaded coin 0.993%, fair coin 0.007%
head 37  loaded coin 0.997%, fair coin 0.003%
head 38  loaded coin 0.999%, fair coin 0.001%
head 39  loaded coin 0.999%, fair coin 0.001%
head 40  loaded coin 1.0%, fair coin 0.0%
head 41  loaded coin 1.0%, fair coin 0.0%
head 42  loaded coin 1.0%, fair coin 0.0%
head 43  loaded coin 1.0%, fair coin 0.0%
head 44  loaded coin 1.0%, fair coin 0.0%
head 45  loaded coin 1.0%, fair coin 0.0%
head 46  loaded coin 1.0%, fair coin 0.0%
head 47  loaded coin 1.0%, fair coin 0.0%
head 48  loaded coin 1.0%, fair coin 0.0%
head 49  loaded coin 1.0%, fair coin 0.0%

that's all

Recommended Posts

Calculate the probability of being a squid coin with Bayes' theorem [python]
Calculate the shortest route of a graph with Dijkstra's algorithm and Python
Calculate the total number of combinations with python
Calculate the probability of outliers on a boxplot
Calculate the regression coefficient of simple regression analysis with python
Calculate the product of matrices with a character expression?
I wrote a doctest in "I tried to simulate the probability of a bingo game with Python"
Save the result of the life game as a gif with python
[Statistics] Grasp the image of the central limit theorem with a graph
[python, ruby] fetch the contents of a web page with selenium-webdriver
The story of making a standard driver for db with python.
The idea of feeding the config file with a python file instead of yaml
Tips: [Python] Calculate the average value of the specified area with bedgraph
The story of making a module that skips mail with python
Create a compatibility judgment program with the random module of python.
Check the existence of the file with python
Search the maze with the python A * algorithm
[python] [meta] Is the type of python a type?
The story of blackjack A processing (python)
The story of making a university 100 yen breakfast LINE bot with Python
[AtCoder explanation] Control the A, B, C problems of ABC182 with Python!
Get the number of searches with a regular expression. SeleniumBasic VBA Python
[AtCoder explanation] Control the A, B, C problems of ABC186 with Python!
[Introduction to Python] How to sort the contents of a list efficiently with list sort
[AtCoder explanation] Control the A, B, C problems of ABC185 with Python!
Hit a method of a class instance with the Python Bottle Web API
Receive a list of the results of parallel processing in Python with starmap
[AtCoder explanation] Control the A, B, C problems of ABC187 with Python!
[AtCoder explanation] Control the A, B, C problems of ABC184 with Python!
How to calculate the volatility of a brand
[AtCoder] Solve A problem of ABC101 ~ 169 with Python
[Python] Get the files in a folder with Python
Prepare the execution environment of Python3 with Docker
[Note] Export the html of the site with python.
[AtCoder explanation] Control the A, B, (C), D problems of ABC165 with Python!
[AtCoder explanation] Control the A, B, C, D problems of ABC183 with Python!
Use logger with Python for the time being
Around the authentication of PyDrive2, a package that operates Google Drive with Python
Make a copy of the list in Python
Check the date of the flag duty with Python
A note about the python version of python virtualenv
[Python] A rough understanding of the logging module
Output in the form of a python array
Convert the character code of the file with Python3
[AtCoder explanation] Control the A, B, C, D problems of ABC181 with Python!
A discussion of the strengths and weaknesses of Python
Get the stock price of a Japanese company with Python and make a graph
[Python] Determine the type of iris with SVM
How to get a list of files in the same directory with python
[Introduction to Python] How to get the index of data with a for statement
How to identify the element with the smallest number of characters in a Python list?
[Python & SQLite] I analyzed the expected value of a race with horses with a win of 1x ②
A memo of misunderstanding when trying to load the entire self-made module with Python3
Python: Calculate the uniform flow depth of a rectangular cross section using the Brent method
I thought about why Python self is necessary with the feeling of a Python interpreter
Extract the table of image files with OneDrive & Python
[Python] A program that counts the number of valleys
A memo connected to HiveServer2 of EMR with python
Calculate volume from the two-dimensional structure of a compound
Learn Nim with Python (from the beginning of the year).
Recommendation of building a portable Python environment with conda