Linear statistical model with two-class variables Suppose I have a set of (discrete) variables, say $X_1,\dots,X_n$. Each $i$ belongs to either class A or class B. When it belongs to A the contribution is $Y_{A,i}f(X_i)$, and when it belongs to B the contribution is $Y_{B,i}g(X_i)$, so that the final outcome is $Z=\sum_{i=1}^N Y_{A,i}f(X_i) I_{\{i \in A\}} + Y_{B,i}g(X_i) I_{\{i \in B\}}$, where $I_{\{.\}}$ denotes the indicator-function.  
The values for $X_1,\dots,X_n$ are given, as are the functions $f$ and $g$. Now suppose I have $M$ obervations for $Z$, what method do you advice to
(i)  Identify whether $i$ is in $A$ or $B$
(ii) Estimate the values for $Y_{A,i}$ (if $i\in A)$ and $Y_{B,i}$ (if $i \in B$).
I tried to estimate the coefficients $Y_{A,i}$ and $Y_{B,i}$ using standard linear regression, assuming that $i\in A$ for all $i$ or $i \in B$ for all $i$, but that gives bad results.
(iii) What if I have $M$ different inputs $X_{j,1},\dots,X_{j,N}$, for $j=1,\dots,M$ and one observation for $Z$ for each input?
Edit:
As an example, take the following R-code:
N=50
M=10
X=rbinom(N*M,1,0.5)-rbinom(N*M,1,0.5) #Generate -1, 0, 1 variables
dim(X) <- c(N,M)
pA=1/2 #pB=1-pA
XA = X #f(x)=x
XB = -X^2 + 1 #g(x)=-x^2+1
isA = as.logical(rbinom(M,1,pA))
YA = rnorm(M,1,5) #generate coefficients
YB = rnorm(M,1,5) #generate coefficients
Z = XA[,isA] %*% YA[isA] + XB[,!isA] %*% YB[!isA] #Calculate Z without noise
Z = Z+rnorm(N,sd=sqrt(var(Z))) #add noise 

Thus given the observations Z and the input data XA, XB (these are known) I want to estimate isA, YA[isA] and YB[!isA].
 A: As you stated the problem, the variables $X_i$ are ordered, so if an $X_i$ is associated with class $A$ for one input set, $X_i$ is also associated to the same class in all the other input sets. This helps enormously to restrict the number of possible models. So your model is of the form:
$$Z=\sum_{i=1}^N Y_i\,f(X_i) I_{\{i \in A\}} + Y_{i}\,g(X_i) I_{\{i \in B\}}$$
Thus, you have a model with $2N$ parameters; $N$ multiplication factors $Y_i$ and $N$ binary variables for the class membership. You can use maximum-likelihood-estimation (MLE) to get the classification and the coefficients. Since you assume that the true $Z$ value is diluted by Gaussian noise, your log-likelihood is (after shifting by an irrelevant constant regarding maximization) proportional to:
$$\text{LL} \propto -\frac{1}{M}\sum_{i=j}^{M}(Z_j-\tilde{Z}_j)^2$$
where $\tilde{Z}_j$ is the model prediction for the given set of observations $X_{j,i}$ and $Z_j$the observed effect for input set $j$ using a specific set of model parameters. Now, the idea is simple: choose the model parameters in such a way that the $\text{LL}$-value is maximized. In the case of 10 predictor variables this can be done by exhaustive search: for every possible class assignment of the $X_i$, solve the linear regression problem. In this way you get for $N=10$ a number of 1024 different models. Pick the model with the highest log-likelihood value. I wrote a small R function that works together with your example to calculate the log-likelihoods:
judgeModels <- function()
{
  resfrm <- data.frame(replicate(M,logical(2^M)),LL=Inf)

  for (i in seq(0,2^M-1))
  { 
    # generate next classification possibility
    isA <- rep(FALSE,M)
    tmp <- i
    for (j in seq(M))
    {
      isA[j] <- tmp %% 2 == 1
      tmp <- trunc(tmp / 2)
    }
    # set design matrix according to classification 
    curXA <- XA; curXA[,!isA] <- 0
    curXB <- XB; curXB[,isA] <- 0
    XAB <- curXA + curXB
    # solve the linear regression problem
    b <- solve(t(XAB)%*%XAB)%*%t(XAB)%*%Z
    # calculate the negative log-likelihood
    resid <- XAB%*%b - Z
    estvar <- 1/length(Z)*sum(resid^2)
    # store the result
    resfrm[i+1,1:M] <- isA
    resfrm[i+1,M+1] <- estvar
  } 
  # order models according to their performance
  # best first, worst last
  resfrm <- resfrm[order(resfrm[,M+1]),]
  rownames(resfrm) <- seq(nrow(resfrm))
  return(resfrm)
}

It returns all possible models together with their log-likelihood score. The lower the score, the better rated the model. Here an example of how to use the function:
resfrm <- judgeModels()
# display top rated models
resfrm[1:5,]
# rank of the true model 
trueidx <- which(apply(resfrm[,1:M],1,function(x) all(x==isA)))
resfrm[trueidx,]

If you play a bit with this function, you will realise that 50 input sets are to less to find the right classification with high certainty. You need a sample size > 1000 to achieve this goal. However, the smaller the variance of the Gaussian noise added to $Z$, the smaller the required sample size to get the right classification, e.g. try Z = Z+rnorm(N,sd=1) and look what happens.
