# What would be sensible weights when fitting simulated probabilities?

I have simulated a set of probabilities $\{p_n\}$ by the hit-or-miss Monte Carlo method. Specifically, my program checks whether or not an $n\times n$ matrix has a certain feature or not; say that, out of $N$ trials, $k_n$ have this feature (success). Then $p_n=k_n/N.$

Question: Now, I want to find the uncertainty (or standard error) $\sigma_{p_n}$ on $p_n$. I'd like to do this so that I can weigh my probabilities when fitting these probabilities as a function of $n.$ My weights are defined as $w_n=\sigma_{p_n}^{-2}.$ This definition seems to be completely standard, but it doesn't go well with my intuition, as explained below. I'd like to know if my intuition is simply wrong (and why), or if there are more suitable weights I could be using.

My problem: As far as I understand, $k_n$ is distributed binomially, so $$\sigma_{p_n}=N^{-1/2}\sqrt{p_n(1-p_n)}$$

This means that a $p_n,$ smaller than some $p_m\leq1/2,$ has a smaller uncertainty than $p_m.$ This makes sense to me (a random walk will not go very far from its starting point if it has a very low probability of taking a step). However, it seems very counter-intuitive to me to weigh a very low probability $p_n$ (where my program perhaps only had $10^2$ successes out of $10^9$ trials) higher than a high probability $p_m$ (which had success, say, half of the time). Using a "relative fluctuation" away from some $p_n$, defined as $$\sigma_{p_n}/p_n=N^{-1/2}\sqrt{\frac{1}{p_n}-1},$$ as my effective uncertainty seems closer to what my intuition tells me is important when choosing my weights.

I'm having a hard time formulating why this is, but it is like this: Flip a coin $100$ times: It comes up heads $54$ times. Flip another coin $100$ times: It comes up heads a single time. Now estimate the bias of each coin. Which estimate are you most confident in? I would say that I'm much more confident in my estimate for the first coin. Should I be? This is the essence of my question.

• Assuming your matrix simulations are independent then $k_n$ does indeed have binomial distribution. However, the success probability ('$p$') is the probability that a realisation of your matrix has the particular feature. Your estimate $\sigma_{p_n}$ would itself be subject to uncertainty, depending on how rare your event is and how big $N$ is etc. Perhaps you can use some elementary probability bounds to estimate the variance of your $\sigma_{p_n}$, – P.Windridge May 17 '17 at 13:02
• Btw, yes it's true that the variance for a binomial with low or high success probability $p$ is smaller than that with moderate $p$. .... consider the limiting cases $p = 0$ or $p=1$. In your example, what is the probability that the true '$p$' for coin B is bigger than $0.1$? Perhaps you're more confident of your estimate of the bias for coin A because you're in a Bayesian world and have a prior belief distribution on the bias :D – P.Windridge May 17 '17 at 13:12
• @P.Windridge 1st comment: I'm assuming that my found $p_n$ are my best guess at a mean, and I therefore use it as the true value in my formulas for $\sigma_{p_n}$. Is this approach flawed? 2nd comment: Perhaps :P, but my point was that coin $A$ has more successes, which just happens to happen around $p=1/2.$ – Bobson Dugnutt May 17 '17 at 13:17
• To help your intuition: with $54$ successes out of $100$ flips, a $95\%$ confidence interval for the underlying probability is (approximately) $[0.45,0.64]$. With $1$ success out of $100$ flips, the interval is (again approximately) $[0.0025, 0.055]$. The uncertainty in the latter is clearly much less than in the former. The relative uncertainty in the latter is much greater than in the former. (That concept is best expressed on an odds scale rather than a probability scale.) – whuber May 18 '17 at 20:43
• You could use either (or even something else). It depends on whether you want to measure the goodness of fit in terms of relative or absolute differences. If it's the former, which might be wise, you should favor a generalized linear model with a Binomial or Poisson response. – whuber May 18 '17 at 21:23

I performed a simple experiment in R where I performed a Monte Carlo simulation to estimate probabilities in the range 0.01 to 0.5 using 50 samples. I then fit the results using two different weighting methods. The first is the inverse of the variance and the second is the inverse of the coefficient of variation. For the second model, I log transformed the data before fitting.

The first model has lower RMSE and the second model has lower MAPE (mean absolute percentage error). If you are trying to minimize the error in your predictions then you may wish to use the first method. If you are trying to minimize the relative error then you may wish to use the second method.

Side note: You may wish to consider an adaptive sample size in your MCS since it is relatively easy to add the estimate of the error and increase the sample size until it is below a desired threshold.

library(tidyverse)

# Monte Carlo Simulations -------------------------------------------------

N = 50
p <- seq(0.01, 0.5, 0.01)
phat <- NULL
for (i in 1:length(p)) {
phat[i] <- sum(rbinom(N, 1, p[i])) / N
}

shat <- sqrt(phat * (1 - phat)) / N

myData <-
data_frame(
p = p,
phat = phat,
shat = shat,
cv = shat / phat,
lb = phat - 1.96 * shat,
ub = phat + 1.96 * shat,
w1 = 1 / shat^2,
w2 = 1 / cv
)

# Plot results ------------------------------------------------------------

ggplot(data = myData, aes(
x = p,
y = phat,
ymin = lb,
ymax = ub,
)) +
geom_point(aes(size = w1)) +
geom_errorbar() +
labs(title = "weighting by inverse of variance",
x = "true probability",
y = "estimated probability") +
geom_line(aes(x = p, y = p), color = "black", linetype = "dashed")

ggplot(data = myData, aes(
x = log10(p),
y = log10(phat),
ymin = log10(lb),
ymax = log10(ub)
)) +
geom_point(aes(size = w2)) +
geom_errorbar() +
labs(title = "weighting by inverse of coefficient of variation",
x = "log true probability",
y = "log estimated probability") +
geom_line(aes(x = log10(p), y = log10(p)), color = "black", linetype = "dashed")

# Weighted least squares --------------------------------------------------

fm1 <- lm(phat ~ p, data = myData %>% filter(shat != 0), weights = w1)

fm2 <- lm(log10(phat) ~ log10(p), data = myData %>% filter(shat != 0), weights = w2)

myData <- myData %>%
filter(shat != 0) %>%
mutate(pred1 = predict(fm1),
pred2 = 10^predict(fm2)) %>%
gather(key, value, pred1, pred2)

myData %>%
group_by(key) %>%
summarize(RMSE = sqrt(sum(value - p)^2/n()),
MAPE = (100/n())*sum(abs((p - value)/p)))

ggplot(data = myData, aes(x = log10(p), y = log10(value), color = key)) +
geom_line() +
geom_line(aes(x = log10(p), y = log10(p)), color = "black", linetype = "dashed") +
labs(title = "ture value vs predictions",
x = "log true probability",
y = "log estimated probability")

• A relevant comparison would be to the unweighted regression. – Hernan May 29 '17 at 7:52