0

I know how to do a standard binomial distribution in python where probabilities of each trial is the same. My question is what to do if the trial probabilities change each time. I'm drafting up an algorithm based on the paper below but thought I should check on here to see whether there's already a standard way to do it.

http://www.tandfonline.com/doi/abs/10.1080/00949658208810534#.UeVnWT6gk6w

Thanks in advance,

James

dotdotdog
  • 1
  • 4
  • Do you know the trial probabilites beforehand? –  Jul 17 '13 at 09:38
  • Yup. And I also know the outcome. – dotdotdog Jul 17 '13 at 10:01
  • Do you mean you want to draw samples from a binomial distribution based on, let's say Trial(1) -> P(Heads)=0.7 and P(Tails)=0.3 and the P(H) and P(T) change for each trial? – Zhubarb Jul 17 '13 at 10:12
  • Not sure if want to 'draw samples'. Using your example, I'd do a bias coin toss n times, with p(H) and P(T) changing on each trial. i would then like to work out the probability of tossing an arbitrary number of heads. – dotdotdog Jul 17 '13 at 10:20

3 Answers3

0

Is this kind of what you are looking for?

import numpy as np

def random_MN_draw(n, probs): # n=2 since binomial
    """ get X random draws from the multinomial distribution whose probability is given by 'probs' """
    mn_draw = np.random.multinomial(n,probs) # do 1 multinomial experiment with the given probs with probs= [0.5,0.5], this is a fair coin-flip
    return mn_draw 

def simulate(sim_probabilities):
    len_sim = len(sim_probabilities)
    simulated_flips = np.zeros(2,len_sim)
    for i in range(0,len_sim)
         simulated_flips(:,i) = random_MN_draw(2, sim_probabilities(i))

    # Here, at the end of the simulation, you can count the number of heads
    # in 'simulated_flips' to get your MLE's on P(H) and P(T).
Zhubarb
  • 11,432
  • 18
  • 75
  • 114
0

Suppose you want to do 9 coin tosses, and P(H) on each flip is 0.1 .. 0.9, respectively. !0% chance of a head on first flip, 90% on last.

For E(H), the expected number of heads, you can just sum the 9 individual expectations.

For a distribution, you could enumerate the ordered possible outcomes (itertools.combinations_with_replacement(["H", "T"], 9))

(HHH HHH HHH) (HHH HHH HHT) ... (TTT TTT TTT)

and calculate a probability for the ordered outcome in a straightforward manner.

for each ordered outcome, increment a defaultdict(float) indexed by the number of heads with the calculated p.

When done, compute the sum of the dictionary values, then divide every value in the dictionary by that sum.

You'll have 10 values that correspond to the chances of observing 0 .. 9 heads.

Gerry

Gerry
  • 1,303
  • 1
  • 10
  • 16
0

Well, the question is old and I can't answer it since I don't know pythons math libraries well enough.

Howewer, it might be helpful to other readers to know that this distribution often runs under the name

  • Poisson Binomial Distribution
imix
  • 1,231
  • 11
  • 13