在R中利用样本量可变的数据拟合二项分布，求解成功概率p

Got it, let's work through this together! You’ve got a collection of binomial trials where each has a different sample size (x = number of attempts) and corresponding successes (y), and you want to estimate the shared success probability p. Here’s how to do this in R, since your data is formatted for that language:

Step 1: Confirm your data setup #

First, let's make sure your vectors are defined correctly (you already have this, but it's good to start here):

x <- c(3, 6, 1, 31, 1, 18, 73, 29, 2, 1)  # Sample sizes (number of trials per group)
y <- c(1, 1, 0, 8, 0, 0, 8, 1, 0, 0)      # Corresponding successes per group

Step 2: Estimate the success probability p #

Since all these trials share the same underlying success probability p, we have two straightforward ways to calculate its maximum likelihood estimate (MLE):

Option 1: Direct MLE calculation (simplest method) #

The MLE for p in this scenario is just the total number of successes divided by the total number of trials across all groups. This comes directly from maximizing the combined binomial likelihood for all your trials.

Run this code to compute it:

total_successes <- sum(y)
total_trials <- sum(x)
p_mle <- total_successes / total_trials
p_mle

When you run this, you’ll get p_mle = 19/165 ≈ 0.115—that’s your estimated success probability.

Option 2: Using a GLM (for future flexibility) #

If you might want to add covariates later (e.g., if p could vary with other factors), using a generalized linear model (GLM) with a binomial family is a great approach. For this case, we’ll fit a model with no predictors (only an intercept) to estimate the overall p:

# Create a response matrix: columns = [successes, failures]
response_matrix <- cbind(y, x - y)

# Fit the binomial GLM with just an intercept
binomial_model <- glm(response_matrix ~ 1, family = binomial)

# Convert the intercept (log-odds) to the probability `p`
p_glm <- plogis(coef(binomial_model))
p_glm

This will give you the exact same result as the direct MLE—plogis() converts the log-odds output from the GLM to a probability.

Step 3: Verify the model (optional) #

If you want to check the model details (like standard errors for your estimate), use the summary() function:

summary(binomial_model)

The intercept’s coefficient is the log-odds of success, and the summary includes standard errors for this value. If you need a standard error for p itself, you can calculate it using the delta method, but for most basic use cases, the point estimate is sufficient.

Why this works #

When you have multiple binomial trials with varying sample sizes but a common success probability, the combined likelihood is the product of each individual trial’s likelihood. The MLE of p ends up being the overall proportion of successes across all trials—this is intuitive and statistically sound for this scenario.

内容的提问来源于stack exchange，提问作者dtcitron