Sorry for the delay, I was going to code this up as an example, but basically my method was going to agree with @Aniko's comment about using a linear regression.
So basically you can run a multiple linear regression where the outcome variable is number of heads and the predictors (or covariates) to the model is the number of tosses of each coin. Doing it so lets you model the scenario with as many coins as you want. Now, more formally, we want to model the following:
$$N_i=Y_{1i}p_1+Y_{2i}p_2+...Y_{Xi}p_X$$
where $N_i$ is the total number of heads on day $i$, $Y_{ji}$ is the $j^{th}$ coin thrown on day $i$, $X$ is the total number of coins, and $p_j$ is the probability of heads associated with the $j^{th}$ coin. So, we would like to be able to obtain estimates of the $p_j$'s.
One way (which is the way we shall illustrate) is to fit the above regression model using least squares. Then, the least square estimators will give the estimated probability of head for each coin (note we exclude the intercept).
To see this method in action, we will simulate data in R
for 5,000 consecutive days with four coins each being flipped a random amount of times each day. Here is the R
code for achieving this:
#The probabilities of heads for the 4 coins
p1 = 1
p2 = .3
p3 = .6
p4 = .45
#The number of days to flip the coins
days = 5000
#Variable for tracking number fo flips of each coin
coin1.flips = rep(NA,days)
coin2.flips = rep(NA,days)
coin3.flips = rep(NA,days)
coin4.flips = rep(NA,days)
#Variable for tracking total number of heads
heads = rep(NA,days)
#Variable for tracking total number of flips
flips = rep(NA,days)
for(i in 1:days){
#Sample a random number of flips for each coin
coin1.flips[i] = sample(10:100,1)
coin2.flips[i] = sample(10:100,1)
coin3.flips[i] = sample(10:100,1)
coin4.flips[i] = sample(10:100,1)
#Flip the coins
coin1 = rbinom(coin1.flips[i],1,p1)
coin2 = rbinom(coin2.flips[i],1,p2)
coin3 = rbinom(coin3.flips[i],1,p3)
coin4 = rbinom(coin4.flips[i],1,p4)
#Calculate total number of heads and flips
heads[i] = sum(coin1)+sum(coin2)+sum(coin3)+sum(coin4)
flips[i] = coin1.flips[i] + coin2.flips[i] + coin3.flips[i] + coin4.flips[i]
}
#Fit the least sqaures model
model = lm(heads~coin1.flips+coin2.flips+coin3.flips+coin4.flips-1)
model
So we set the initial probabilities for the 4 coins to be $p_1 =1$, $p_2=0.3$, $p_3=0.6$, and $p_4=0.45$
Running the code we obtain the following estimates for the four probabilities:
> model
Call:
lm(formula = heads ~ coin1.flips + coin2.flips + coin3.flips +
coin4.flips - 1)
Coefficients:
coin1.flips coin2.flips coin3.flips coin4.flips
1.0025 0.3022 0.6001 0.4459
And so we see that the coefficients (which are the estimated probabilities) are extremely close to the true solutions.