How can I identify probability mass functions with known proportions of contribution in a series of histograms? Suppose I have two real-valued PMFs $f(i)$ and $g(i)$ with $1\le i\le N$ and $N$ is about 1 billion.  The functions $f$ and $g$ can be assumed to be continuous on the positive integers, one-sidedly so at the endpoints of course. (Meaning if i and j are close, then $f(i)$ and $f(j)$ are close.)
Lets say I have several samples of more than 10 billion observations each.  I tally each sample to get a series of histograms. The histograms basically give us functions $s_1(i)$, $s_2(i)$, $s_3(i)$ and so on.
For each sample, I know the percentage of the observations that come from $f$ and the percentage that come from the distribution described by $g$.
Can I use the samples to reconstruct the distributions $f$ and $g$.
My first guess was to do linear regression.  Although, I'm unsure if linear regression makes use of all the information.
 A: The challenge of this question lies in the large amount of data.  Let's approach a solution by making a succession of approximations, each faster than the last, until a solution becomes practicable.  That solution indeed is a set of regression problems, as indicated in the question.  A method to speed it up is proposed and implemented, resulting in a four order-of-magnitude speed increase.
The results obtained with simulated (yet realistic) data should be sobering: they point to inherent limitations in the accuracy that can be expected when estimating $2N$ parameters from approximately $2N\times 10\times k$ observations of extremely rare outcomes.
Theoretically ideal solutions
Use Maximum Likelihood or Bayesian procedures to estimate the $2N$ parameters $f_1, \ldots, f_N$ and $g_1, \ldots, g_N$.  But we don't want to have to optimize a likelihood that depends on $2N \approx$ two billion linked variables!
An excellent approximation
Sample $j$, $j=1, 2, \ldots, k$, is a realization $(x_1^{(j)}, \ldots, x_N^{(j)})$ of a random variable $(X_1^{(j)}, X_2^{(j)}, \ldots, X_N^{(j)})$ which is a mixture of two multinomial distributions $F$ (with probabilities $f_i$) and $G$ (with probabilities $g_i$) in known proportions: let $\alpha_j$ be the proportion of $F$,whence $1-\alpha_j$ is the proportion of $G$.   With $N$ so large, the $X_i$ are for all practical purposes independent and they have Poisson distributions.  This splits the problem into a billion subproblems, one for each index $i$:
Given that $X_i^{(j)}$ has Poisson distribution with expectation $\alpha_j f_i + (1-\alpha_j) g_i$, $j=1, \ldots, k$, estimate the two parameters $f_i$ and $g_i$.
In more familiar notation this says
$$X_i^{(j)} \sim\text{Poisson}(\mu_i^{(j)}),\quad \mu_i^{(j)} = \beta_0 + \beta_1 \alpha_j$$
with $\beta_0 = g_i$ and $\beta_1 = f_i - g_i$.  This is a Generalized Linear Model (GLM) for $X_i$ with a Poisson distribution and identity link.  To make the estimates into an actual probability distribution, we can clamp them to the interval $[0,1]$ and rescale them all to make the $\hat{f_i}$ and the $\hat{g_i}$ separately sum to unity.
The problem with this is that running a billion GLM estimates would take a full week of computing (in R).
Speeding things up
We can approximate a Poisson regression by means of a Weighted Least Squares regression: because the variance of a Poisson$(\mu)$ variate equals $\mu$, we can use the observed values $x_i^{(j)}$ as weights.
Unfortunately, because the weights vary with $i$, this still requires conducting a billion separate linear regressions. An R implementation goes only twice as fast as the GLM solution.
Another approximation
If we drop the weighting, the solution is not quite as good as before, but it works pretty well nevertheless.  Without weights, notice that every one of the regressions uses the same explanatory variable $\alpha$. Because the Ordinary Least Squares estimator is linear, we can do this regression once and for all and then each estimate is just a simple, fast, linear combination of the data.  This is extremely fast: whereas the original GLM solution requires 3 to 4 seconds for $N=1000$, the same machine can perform $2\times 10^7$ OLS regressions in the same time.  Extrapolating this performance suggests cases where $N$ is a billion and $k$ is small can be computed in less than ten minutes (using a single core on a PC: doing this on a GPU could easily drop the time to less than one second plus I/O overhead).
Code
This working R code illustrates the ideas and compares the GLM to the OLS results.  First we need to generate some potentially large datasets with multinomial distributions, after stipulating values for $N$, $k$, the totals in the datasets, and $\alpha_j$:
set.seed(17)                        # Makes output reproducible
n <- 10^3                           # Number of indexes to f and g
k <- 3                              # Number of samples
alpha = seq(0.3, 0.7, length.out=k) # Sample proportions (of f)
totals <- rep(10 * n, k)            # Sample totals, per sample

# Parameter and data generation
f <- rgamma(n,1,1); f <- f / sum(f)
g <- rgamma(n,1,1); g <- g / sum(g)
x <- sapply(1:k, function(i) rmultinom(1, totals[i], alpha[i]*f + (1-alpha[i])*g))

Here is R's implementation of the GLM solution.
system.time({
  raw <- apply(x, 1, function(y) {
    coef <- glm(y ~ alpha, family=poisson(link="identity"), start=c(1/n, 1/n))$coeff; 
    cumsum(coef) # Assemble hat(beta_0) and hat(beta_1) into hat(f) and hat(g)
  });
  z <- apply(raw, 1, function(y) {x <- pmax(0, y); x / sum(x)}) # Normalize the fits
})


  user  system elapsed 
  3.44    0.03    3.47 

There were 50 or more warnings (use warnings() to see the first 50)

The warnings are expected when using a linear link for Poisson regression; as we will see, they do little harm.
The OLS solution requires precomputation of the regression solution, implemented by lm.init:
lm.init <- function(alpha, ...) {
  k <- length(alpha)
  i = matrix(c(rep(c(1, rep(0,k)), k-1), 1), k)
  apply(i, 2, function(y) lm(y ~ alpha, ...)$coeff)
}
system.time({
  raw.ols <- x %*% (t(lm.init(alpha)) %*% matrix(c(1,0,1,1), 2))
  z.ols <- apply(raw.ols, 2, function(y) {x <- pmax(0, y); x / sum(x)})
})

The timing for $N=1000$ is $0$ seconds; for $N=10^7$ it is $1.6$ seconds.
Let's compare the estimates to the true values for the GLM and OLS cases:
par(mfrow=c(2,2))
plot(g, z[,1], ylab="g.hat", main="GLM"); lines(c(0,max(g)), c(0,max(g)), col="Red")
plot(f, z[,2], ylab="f.hat", main="GLM"); lines(c(0,max(f)), c(0,max(f)), col="Red")
plot(g, z.ols[,1], ylab="g.hat", main="OLS"); lines(c(0,max(g)), c(0,max(g)), col="Red")
plot(f, z.ols[,2], ylab="f.hat", main="OLS"); lines(c(0,max(f)), c(0,max(f)), col="Red")


Although ideally the estimates would track along the diagonal lines of equality (indicating the points where no error is made), this amount of scatter is unsurprising: in the simulation the counts averaged only $10$ with a variance around $10$, the proportions $\alpha_j$ did not change much (they ranged only from $0.3$ to $0.7$, and only $k=3$ samples were observed.  The symmetric scatter around the lines of equality demonstrates that these estimates are unbiased.
As another check, we may compare the GLM and OLS results.  These are shown on square-root scales to resolve what happens with the tiniest parameters.  The GLM estimates are plotted on the horizontal axes and the OLS estimates on the vertical axes:
plot(sqrt(z[,1]), sqrt(z.ols[,1]), main="g")
plot(sqrt(z[,2]), sqrt(z.ols[,2]), main="f")


The two methods are actually pretty close: OLS is not a bad approximation to Poisson regression. The approximation is worst for the small probabilities, as expected.
